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FOREWORD 

The preparation of this report has been sponsored by the National 
Aeronautics and Space Administration under Contract No, NASr-109. The 
initial development of this computer program was sponsored under Contract 
No. AF 49(6 3 8) - 792 by the United States Air Force Office of Scientific Research, 
under Contract No. AF 40(600)-804 by the USAF, as part of Project Defender 
by the Advanced Research Projects Agency, Department of Defense, and by 
Cornell Aeronautical Laboratory Internal Research, Modifications of the 
original program have been made in connection with research performed 
under Contract No. NASr-109 for the NASA and Contract No. AF 33(657)-8860 
for the Aerospace Research Laboratories of the USAF, 

The original version of the program was coded for use on an IBM 704 
computer by Dr, John Fleck, Department Head, and Duane Larson of the 
Computer Services Department of CAL. The program was recoded in 
FORTRAN IV language for use on an IBM 7044 computer by Mrs. Camille 
Fiore, also of the Computer Services Department. 
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ABSTRACT 


Q <o1i 6 

This report describes a computer program developed at CAL for the 
numerical solution of quasi-one -dimensional, inviscid expansions of reacting 
mixutres, The analytical techniques employed and results obtained using 
the program have been previously reported. The details of the coding and 
operation of the program are discussed in this report. Program cards which 
have been written in FORTRAN IV language for use on an IBM 7044 computer 
are available upon request. 

The computer program in its present form is capable of handling a 
general gas mixture, including up to 20 chemical species and 64 reactions. 

The species are assumed to undergo vibrational and electronic excitation in 
addition to coupled chemical reactions. The vibrational and electronic degrees 
of freedom are assumed to remain in thermodynamic equilibrium while the 
chemical reactions may be assumed to be either frozen, in equilibrium, or 
to proceed at finite reaction rates. Additional options in the program are 
provided for treating ionized nozzle flows, assuming the degree of vibra- 
tional excitation to be frozen at the reservoir value, or including the effects 
of moderate virial imperfections on high-density nozzle flows. 

While the basic program computes the properties of an expansion from 
an equilibrium reservoir state through a nozzle of specified geometry, modi- 
fications of the program have been made to permit specifying the stream- 
tube pressure variation rather than an area distribution and to permit a 
finite- velocity, equilibrium or nonequilibrium initial state. The program 
cards for these modified versions of the program are also available upon 
request. 
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1. INTRODUCTION 


In high temperature expansions of reacting gas mixtures the flow transit 
times are usually not sufficient to maintain equilibrium. Departures from 
equilibrium in turn affect the gasdynamic properties of the flow. Thus, in 
order to define the flows in high-temperature, hypersonic wind tunnels and 
rocket nozzles, and around hypersonic vehicles, the nonequilibrium effects 
on the fluid -mechanical properties of the flows must be taken into account. 

In the course of studying nonequilibrium phenomena in high-speed flows a 
computer program has been developed at CAL to obtain numerical solutions 
for quasi-one -dimensional, inviscid expansions of chemically reacting 
mixtures. ^ This program, which is also capable of generating solutions for 
the limiting cases of frozen and equilibrium flow, has proved to be a useful 
tool, not only in the above-mentioned situations, but also in the design and 
interpretation of fundamental experiments on nonequilibrium flows. 

The purpose of this report is to describe the methods of solution 
contained in the program and the use of the program. The report is intended 
to accompany a copy of the program cards, which are available upon request. 
The program cards are in FORTRAN IV language and are coded for use on an 
IBM 7044 computer. This report is similar in purpose to Reference 2 which 
describes the program developed at CAL for computing the nonequilibrium 
flow behind normal and bow shock waves. 

Several authors have reported numerical solutions for nonequilibrium 
expanding flows (e.g. Refs. 3-6). The initial work in this area at CAL was 

7 

performed by Hall and Russo who obtained exact solutions for the nonequilib- 
rium nozzle flow of a diatomic gas and a diatomic gas in an inert diluent. 

g 

Then Boyer, Eschenroeder, and Russo reported equilibrium and approximate 
nonequilibrium solutions for nozzle expansions of air over a wide range of 
reservoir conditions. The approximate nonequilibrium solutions were obtained 
using the approximate freezing criterion of Reference 7 for the oxygen 
dissociation reaction. Then, the method of exact numerical solution for 
nonequilibrium expansion flows, which is contained in the computer program 
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F urther 


described here, was presented by Eschenroeder, Boyer, and Hall. 

development of this method, particularly for near equilibrium flow, was 

9 

reported by Mates and Lordi. Applications of the computer program 

employed in the method of Reference 1 have included hypersonic wind tunnel 

flows of high- tempe ratur e air, ^ rocket-nozzle flows of hydrogen-carbon 

1 Z 1 3 

and hydrogen-oxygen mixtures, ’ and body streamtube flows of high- 
14- 17 

temperature air. 

The basic nozzle-flow computer program has the capability of obtaining 
numerical solutions for inviscid expansions of a reacting mixture from an 
equilibrium reservoir state through a streamtube of arbitrary geometry. The 
vibrational and electronic degrees of freedom are assumed to remain equilibrated 
with the translational and rotational modes; hence, the effects of vibrational 
and electronic relaxation on the chemistry are not included. The chemistry 
may be assumed to either be frozen, remain in equilibrium, or proceed at a 
finite rate in the expansion. 

The method by which the composition and thermodynamic properties are 
found for a given mixture in chemical equilibrium at a specified temperature 
and pressure is discussed in Section Z. 1. The extension of this method to 
describe an equilibrium expansion from the specified reservoir condition is 
discussed in Section Z. Z. The computation for an expansion with frozen chemical 
composition is described in Section Z.3. 

The procedure employed in the computer program for obtaining non- 
equilibrium expansions from a given reservoir condition is covered in Section 
3. In Section 3.1, the general approach followed in solving the governing 
equations is described. The problem is basically one of obtaining the numerical 
solution for a system of coupled, first-order nonlinear differential equations. 

Some difficulty is encountered in starting the solution in the reservoir 
because the gradients of the flow properties and composition are assumed 
to vanish there. The method employed to start the solution is to consider the 
nonequilibrium flow in the initial portion of the expansion to be a perturbation 
about the equilibrium solution. The analytical technique for locating a point in 
the expansion for starting the numerical integration is described in Section 3. 2. 



The modification of the technique of solution when the numerical integration 
is started upstream of the nozzle throat is covered in Section 3. 3. The 
perturbation calculation used in the starting procedure is described in Section 
3.4, Finally the numerical integration scheme used in the program and some 
of the numerical difficulties encountered in seeking nonequilibrium solutions 
are discussed in Section 3. 5. 

The method for specifying the kinetic model to be used in a given 
calculation is presented in Section 4. It is necessary to provide thermo- 
dynamic data for the specific enthalpy and chemical potential of each species 
included in the mixture. Also, the forward rate constants must be given for 
each reaction included in the gas model. The program which accompanies 
this report has a capacity for a gas model comprised of ZO species and 64 
reactions. This capacity is determined only by the storage limitation. 

Certain options in the program could be eliminated if larger models need be 
considered , 

The following alternative options, discussed in Section 5, are also 
included in the program. Nonequilibrium nozzle flows of ionized gases are 
considered by treating electrons as a separate chemical element. The 
vibrational mode can be assumed frozen at the reservoir degree of excitation 
while the chemical reactions are assumed to be frozen, in equilibrium or in 
nonequilibrium. Also, moderate gas imperfections may be accounted for in 
the initial near equilibrium portion of the expansion. 

Instructions on the preparation of input data for the computer program 
are given in Section 6. Also, the calculation of a nozzle expansion of air is 
used to demonstrate the use of the program. Some of the results of this 
calculation are used to describe the form of the output data and to relate some 
of the experience gained in the use of the program. The results of the sample 
calculation also serve as a test case for the program. 

In Section 7 a description of the function performed by each subroutine 
in the program is given. Flow charts are provided for the more complicated 
subroutines. These descriptions of the subroutines are intended as a guide 
in relating the analysis and governing equations presented in the report to the 
FORTRAN IV program. Also, a complete list of the FORTRAN symbols used 

in the program including a description of each is given in Section 7. 
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Modifications of the basic program have been made in order to treat 
other applications. One modification permits the solution of a nonequilibrium 
expansion from an equilibrium stagnation point through a specified pressure 
variation. This version of the program is applicable to the flow from a 
stagnation point through the streamtube nearest a body for which the surface 
pressure variation is known. The second modification of the program lends 
itself to the study of quasi -one - dimensional flows starting at a known non - 
equilibrium state and proceeding through either a specified area or pressure 
variation. This modification can be used to compute the relaxing flow behind 
shock waves (constant area streamtube) or along a streamline in the outer 
portion of a blunt-body flow field (specified pressure distribution). Yet 
another variation of the program is aimed at solving the nonequilibrium 
expansion from an equilibrium, but finite- velocity initial state. The 
application of this version of the program is to non- reflect ed -type shock 
tunnels or to the flow downstream of an MHD accelerator. The governing 
equations and the associated changes in the use of the program for these 
three modifications are discussed in the Appendix. The FORTRAN IV language 
program cards for each of the three modified versions of the program are 
also available upon request. 
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2. METHOD FOR COMPUTING EQUILIBRIUM RESERVOIR STATE, 
AND EQUILIBRIUM AND FROZEN FLOWS 


The computer program calculates the equilibrium composition of a 
reservoir state for a specified temperature and pressure. Then the variation 
in composition and flow properties in an equilibrium expansion from the 
specified reservoir state can be computed. Also, the solution can be obtained 
for an expansion in which the composition is frozen at the reservoir value. 

The methods employed in these computations are described in this section. 


2, 1 Reservoir Calculation 

The determination of the composition of a specified gas mixture for a 
given equilibrium state proceeds as follows. If there are c chemical elements 
and chemical species in a multicomponent mixture the chemical formula 
of the a. th species may be represented as 

- fat > °v 2 ’ °Vc) * = l > 5 u) 


The specification of the species to be included in the calculation is then through 

the matrix oc . - , where oc • • is the number of atoms of the J th element per 

*4 * — 

molecule of the 2 th species. 

There are (s-c) linearly independent relations which may be written as 
equilibrium formation reactions for the species in terms of the elements 


4 = c +! > c + 2 , 


> 5 


(la) 


* Z M, 

t Sl 

It is not necessary that the c chemical elements be chosen as independent 
species from which the other species are formed. In fact, the chemical 
elements are not necessarily included as separate species in the mixture. 
However, the independent species or components must be linearly independent 
combinations of the chemical elements. In other words, the rank of the 
oc • ■ matrix must be c . 
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As will be discussed in ensuing sections it is sometimes necessary to 
choose species other than the chemical elements as components. The species 
chosen as components are listed in the first c rows of the o c . matrix. The 
computer program then writes the formation reactions for the ( s-c ) dependent 
species in terms of the c components as 


c 




a = c+ I 9 c + 2>’* • > s ( 2 ) 


In addition to specifying the chemical species and the chemical elements 
to be included in a computation, the number of gram-atoms of each chemical 
element contained in the mixture must be given. The total number of 
gram-moles of the components which are present in the mixture is given by 

c 

This relation now specifies the element composition of the mixture in terms 
of the components rather than the chemical elements. 


As shown in Reference 8, enforcing mass conservation in each of the 
formation reactions and global mass conservation lead to the following rela- 
tion for the mole fractions of the components, X* , in terms of the <z' and 

t fj 

the mole fractions of the other ( 5 _ c ) species 4 


*/ [V^o] x. 


where 



~ C +■ * * » 5 


(4) 
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At chemical equilibrium the dependent species concentrations are related to 
the concentrations of the components through the equilibrium constants for 
the formation reactions (Eq. (Z) ). 


x / - ** r 


fi, -0 


tt x >/ y = c +i , c + z> ■ • •? s (5) 

i-t * 


where 


*>, - ^ 


x 






t 


( 6 ) 


Equation (5) may be used to reduce Equations (4) to a set of c simul- 
taneous equations for the concentrations of the independent species or 
components . 

From the above development the calculation of the reservoir conditions 
for a given temperature and pressure is seen to require specification of the 
species included in the mixture, the chemical potential of these species, and 
the relative number of gram-atoms of each chemical element present. 

Having specified these items the ^.Vare used as initial estimates for the 
independent species concentrations. Using the condition of chemical equilibrium 
(Eq. (5) ) the dependent species concentrations are computed and then the 
results are used to test whether Equation (4) is satisfied. 

A Newton-Raphson iteration scheme is used to solve Equation (4) in the 

form 

*7 (*,’■■■’ * f ) -- 9, - 1 & ' f> ft x* - x / - 0 

( 7 ) 

i - l,2> ■■■ c 

where again the s-c dependent concentrations X ■ > - c +/,5 are found 
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from Equation (5). Newton's method for finding the (si + t ) correction to the 
estimates for X • employs the linear term in a Taylor expansion about the 

* /t 

previous value )(, . Enforcing the condition that F. - 0 for the ( / ) 

J ^ 

iteration, the fractional corrections for the a , & values of the X • are found 

from 



In the program these equations are written in the form 



Letting Jl denote the fractional correction 

AV / A/ 



the (/ i + l ) estimates are obtained from 


A.+ ( 

x. 



^ “ ly * ' * > c 


(9) 


If at a given step in the iteration, the (Aft ) guess for X - found from Eq. (9) i 

a + r J - 

zero or negative, then X. is found from 


>t + 1 



(9a) 


Once the species concentrations for the specified reservoir temperature 
and pressure have been found, the mixture molecular weight at the reservoir 
condition is calculated from 


s 



( 10 ) 
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Then, the reservoir density may be computed from the equation of state 


= p'' (X X *%) 


( 11 ] 


For a mixture of ideal gases the equation is 

-p' = TgX 7 ^ ( 1 1 a) 

r Tri 

Finally the specific reservoir or stagnation enthalpy is computed from 

s 

H =X Jvj (12) 

r 1 

where, for a mixture of ideal gases, the species enthalpies, ^ . , are functions 
only of temperature. In order to complete the reservoir computation both 
the species enthalpy and chemical potential must be specified. The method 
employed in the program for specifying these data is discussed in Section 4. 


Z. Z Equilibrium Expansions 

In the above formulation, the determination of the equilibrium reservoir 
state was seen to require the condition of chemical equilibrium (Eq. (5)), the 
conservation of mass in the formation reactions (Eq. (4)^, and global mass 
conservation. In addition, an equilibrium expansion requires the conservation 
of momentum and energy. For quasi -one -dimensional, inviscid flow the 
global conservation of mass, momentum and energy may be expressed as 


p u. A = y» 

= constant 

(13) 

u du. + dp 

= 0 

(14) 

H + y u 1 = 

= constant 

(15) 
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Now, the entropy of a chemically reacting mixture is 


A- -L- 
R a ‘ % 


X 




5 5 



where 


s 0 _ — i 

/ T 


( 16 ) 


(17) 


The entropy change for a chemically reacting mixture may be written 


18 


( 18 ) 


Then the condition of chemical equilibrium (Eq. (5)), may be rewritten as 

1 8 

follows by enforcing a maximum in the free energy 



(19) 


Also from Equations ( 14) and ( 15) 


<tN 




( 20 ) 


for a one -dimensional, inviscid flow. From Equations (18), (19), and (20) the 
inviscid flow of a mixture in chemical equilibrium is seen to be isentropic. 
Thus, an equilibrium expansion from a given reservoir state can be solved 
in the same manner as the reservoir state by adding the constraint that 
successive points in the expansion constitute an isentropic path. 

The equilibrium solution is obtained by taking successive temperature 
steps from the reservoir value. The results of the previous step are used as 
initial estimates for the pressure and concentrations of the components, and 
the Newton-Raphson method is employed to find the correct pressure and 
composition at the specified temperature and entropy. The governing equation 
which must be satisfied by the C components and the pressure is 
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F . - 0 

1 


^ ~ / y * “ } C > C + \ 


where F- for =■ /,*** ? c is given by Equation 


(7) and 


f., -Z x, (s; - % s 

i = f 


( 21 ] 


Equation (21), which specifies the isentropic path, provides the additional 
equation needed to find the pressure at successive temperature steps. 

As cited above the equilibrium composition and pressure can be computed 
at successive temperature intervals. The other flow properties can then be 
found at each temperature as follows. Once the iteration for the pressure and 
composition converges, the local molecular weight is computed from Equation 
(10). Then the density is obtained from the equation of state (Eq. (11a)). If 
T', , and p ' are nondimensionalized by the reservoir values the state 

equation may be written 

p -2— 

r nq nit) 


The static enthalpy for the mixture is computed from Equation (12) and the 
velocity from the energy equation (Eq. (15)). The value of the Mach number 
at each point is obtained from the following formula. 



( 22 ) 


where ^ and are the values of and p at the previous computational step. 


In the above procedure the solution for an expansion in chemical 
equilibrium is obtained as a function of temperature. In a nozzle-flow 
expansion it is desired to obtain the solution as a function of area ratio. Here 
the following procedure is employed to find the equilibrium-flow solution as 
a function of area ratio. The equilibrium solution is obtained at successive 
temperature intervals, computing the flux,^^, at each point. The maximum 
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pu. and the corresponding temperature determine the throat conditions and 
the critical mass flow for the equilibrium expansion. Having thus found the 
mass flow the area ratio corresponding to each temperature can be found 
from Equation (13). 

In the computer program the throat, i. e, , the maximum pu, is located 
and then the equilibrium-flow calculation is restarted at the reservoir. The 
solution is then obtained at fixed temperature intervals, computing the area 
ratio at each point. Since the temperature intervals become very small when 
see ki n g the critical mass flow, the results of that calculation are not printed 
except for the throat values. Relatively few steps are repeated by this 
restarting of the calculation at the reservoir. 

The maximum value of pu . is found by the following method. At each 
temperature step the value of pu- is compared with the value at the previous 
two steps. Denoting the values of pu as f , > an< 3 i n order of 

decreasing temperature, the computation is continued at the original interval 
in temperature until . Once pu passes through a maximum, the 

temperature step size is refined to locate the maximum precisely. 


The refinement of the temperature step proceeds on the following basis. 

If f - r , the temperature interval is halved and the computation restarted 
3 ~ T j 

at the temperature corresponding to .If ^ t ^ ie computation is 

resumed with the smaller temperature interval starting at the temperature 
corresponding to f . The temperature step is continually reduced whenever 
f ^ £ , When the temperature interval becomes less than a specified value 

the sharpness of the maximum is examined. If fc* *,-**,} !/V, is less 

than .00001 the throat temperature is set equal to that corresponding to the 

current value of f , If the above condition is not satisfied a parabolic fit 

2 

for fCr) is obtained from the current values of , f ^ , f . 


n - 0 - JjLL- (t-T ) + IllhlAi i . (j-T'f 
+ Z(AT) 2 (AT')' { %) 


(23) 


where AT is the size of the temperature step and T, , > and T 3 are the 

temperatures corresponding to ^ , and-f 3 . The throat temperature is 

then taken as that value for which cLf/dT obtained from Eq. (23) vanishes. 
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( 24 ) 


T* -T + ft' filAI 

* 2(7; ff 3 - 2 /*) 

The equilibrium solution is obtained in intervals of T of . 0 1 , The 
solution can be terminated at any specified temperature below the throat 
temperature. Since intervals of .01 are taken, the specified temperature 
stop must be Ss . 0 1 because the Newton-Raphson iteration procedure 
diverges for T 0 > within the digital accuracy* If the equilibrium solution 

were desired at temperatures below . 01 J / then the program could be 
modified to decrease the step size below a given temperature. 


2 . 3 Frozen FI o w 


The computer program also contains the option of generating the 
solution for an expansion in which the chemical composition is frozen at the 
reservoir state. As may be seen from Equation (IB) an inviscid, frozen 
expansion is also isentropic. Using the value obtained for the reservoir 
entropy the frozen flow may be determined by direct algebraic computation. 


The procedure employed for the frozen flow solution is to specify a 
temperature interval below the reservoir temperature. For the given 
temperature and yu.°j are calculated for each species. Then using the 
reservoir entropy and composition the pressure is obtained from Equations 
( 1 6) and ( 17 ). 


r 


A- ^ 


V. 


- A TIL 


X K,.j 

If * 1 


- 


(25) 


Next the density and static enthalpy are calculated from the thermal and 
caloric equations of state (Equations (lib) and (12)). The velocity may then 
be evaluated from the energy equation (15) using the reservoir, or stagnation 
enthalpy. 


By following the above procedure the frozen flow solution is obtained 
at specified temperatures. The solution is found in terms of the nozzle area 
ratio by the same method employed for the equilibrium solution. The 
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maximum in pu. is found in order to locate the nozzle throat. The area ratio 
is then found at each computational step using the continuity equation. Also 
a Mach number is computed at each step using Equation (22). 
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3. NONEQUILIBRIUM EXPANDING FLOWS 


In this section the nonequilibrium flow computation is discussed. First 
the analysis, governing equations, and the general method employed in the 
solution are discussed. Then the method for starting the numerical integration 
and the perturbation calculation used in the starting procedure are described. 
Finally the numerical integration scheme used in the program is presented. 


3. 1 Analysis and Governing Equations 


The fluid considered is a multicomponent mixture composed of S 
chemical species. These species undergo /c coupled, chemical reactions 
of the form 



where >). . and >). . are the stoichiometric coefficients of the reactants and 
products and Jh and Jk K . are the reaction-rate constants for the forward 

‘ A* 

and reverse reactions, respectively. Notice that these reactions differ from 
the formation reactions (Equation (1)) and that here /t is not necessarily 
equal to s-C . However, the equilibrium formation reactions may appear 
in the nonequilibrium reaction mechanism. 


The nonequilibrium flow is governed by the conditions of conservation 
of mass, momentum, and energy (Equations (13)-{15)), the thermal and 
caloric state equations (Equations (11) and (12)), and the species conservation 
equations. Again the enthalpy and chemical potential for each species must 
be specified. For the nonequilibrium flow the condition that the chemical 
elements be conserved is written 

s 


I 

=/ 



*- 




(27) 
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or alternatively 




OC 


i* 


■ */ 




= 0 


JL =1,2, ■ • • c (28) 


In the nonequilibrium calculation the equations are written using the concen- 
tration variable /. which is in units of moles per gram. 

3 

The relation between the concentration variable Xj used in the equili- 

d 

brium flow calculation and that used in the nonequilibrium calculation, 
is 



(29) 


Since the species gradients along the nozzle are connected by the element 
conservation equation (Eq. (28)) there are ( s-c ) independent species con- 
servation equations given by 

A 

j = c +i } c +Z, • • • > s ( 3 0 ) 




d X, 


where 



(31) 



(32) 


(33) 


In the above equations rC is the axial distance along the nozzle, <*: ' , divided 
by the characteristic length £ . Also, R and R are the forward and 
reverse rates of the Ji th reaction. 

Due to the complexity of the chemical production terms on the right 
hand side of the species conservation equations, the nonequilibrium solution 
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is obtained numerically. However, since all the flow gradients are assumed 

zero in the reservoir the integration cannot be started there. In the case of 

nonequilibrium solutions for diatomic gases the integration has been started 

7 

using a series expansion in inverse powers of the nozzle area ratios. This 
procedure would be quite complicated for a multicomponent mixture. How- 
ever, in the initial portion of the nozzle, the solutions for a diatomic gas 
indicate that the flow remains very near equilibrium. Thus, the procedure 
used to start the numerical integration in the present case is to treat the 
initial portion of the nonequilibrium solution as a perturbation about the 
infinite-reaction- rate equilibrium solution. The computation of the pertur- 
bations is discussed in Section 3. 4. The procedure for starting the numerical 
integration and the governing equations for the nonequilibrium computation 
are discussed below. 


The nonequilibrium expansion is governed by 5 + 3 first order differential 
equations, element conservation, Equation (27), species conservation, 

Equation (30), and the conservation of mass, momentum, and energy. For 
the nozzle-flow problem the area distribution is specified and the unknowns 
are the species concentrations, u , T , and p . The thermal and caloric 
state equations are used to relate p and H to these dependent variables. In 
the computer program for the numerical solution of this problem the governing 
equations are arranged as follows. The energy and state equations are used 
to eliminate the pressure and velocity derivatives in the momentum and 
continuity equations. Differentiating Equations (13) and (15) and eliminating 
du fd * yields 

<jlrvA I dH 

d d Cl * d 

Using Equation (12) to rewrite dti/d? yields 






cc 1 d£n,p ii 1 d jtw A 

7r( c d# 


(34) 
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whe re 


dJL; 

c - — 1 ~ 

dT 


( 35 ) 


Now differentiating Equations (11) and (15) and substituting the results in 
(14), the momentum equation becomes 

4H T7n. . a?. <CT. T *2L , 0 

d* * Tq 4* '*1 ** ’’’l"' 

Using Equation (12) and the fact that 




leads to 




i $ 

Using Equation (34) this result becomes 

C ifj j ^- 2 - ^ ^ 7 ” / I # * \ U^rup U dbi A ^ 

4r ^T/ att ~ %T 


(36) 


For a prescribed geometry, Equations (28), (30), (34), and (36) are a set of 

dY; dT i d£*vp 

S + Z equations for the unknown slopes ^ > and ^ at a 


point in the numerical integration. The quantity dJbvp/d * is not eliminated 
from Equations (34) and (36) because of the procedure adopted for starting 
nonequilibrium solutions upstream of the nozzle throat. 


As discussed in more detail in Section 6> the nozzle geometry may be 
specified in either of two ways. For a wedge, cone, or hyperboloid geometry 
the area ratio at a given nozzle station is given by a simple quadratic in ^ . 
Different shapes may be specified for the converging and diverging sections. 
Accordingly, for < 0 the area ratio is given by 


A(d) = A,t + A 3 d* 


(37) 


and for /X > 0 

a (* ) - 6, + 8** + BjX 1 


(37 a) 
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The other method of specifying the nozzle geometry is applicable to 
contoured shapes. In this case the nozzle area-ratio variation may be fitted 
with a polynomial in seven different intervals, two in the converging section 
and five in the diverging section. For reasons which become evident in 
Section 3. 2 the polynomials in the first intervals in both sections must be 
quadratics. The polynomials in the other regions may be up to sixth order. 

Once the state of the flow is known at a given nozzle station d and 

— — is obtained from the specified , slopes dt/d # , dT/dy, » 

and dl*,p / d 4 , maybe calculated from the above set of s+ 2 differential 
equations. These equations may then be integrated numerically to obtain 
the temperature and composition at a new nozzle location. The numerical 
technique employed is discussed in Section 3. 5. Once a new temperature 
and composition have been determined the enthalpy may be calculated from 
Equation (12), the velocity from Equation (15), the density from Equation (13), 
and the pressure from Equation (11), At each point in the nonequilibrium 
solution a Mach number is calculated based on a "speed of sound" given 
by 

Notice that the above defined speed of sound is not in general the speed of 

1 9 

propagation of small disturbances in a relaxing medium. Also, at each 
point a value of the entropy is computed from Equation (16). 

In evaluating the density from the continuity equation (Equation (13)) 
it is tacitly assumed that > the critical mass flow is known. If the flow is 
appreciably out of equilibrium upstream of the nozzle throat, -m is not known 
at the outset of the calculation. The method for handling this situation is 
discussed below together with the description of the procedures for starting 
the numerical solution. 

The nonequilibrium flow solution may be stopped at either a specified 
nozzle station or at the same temperature which is used for the frozen and 
equilibrium solutions. Since the nonequilibrium solution may be time con- 
suming (i, e. longer than 15 minutes on an IBM 7 044) for a complex, chemical- 
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kinetic model, provision can also be made for saving the computation at a 
given point, as described in Section 7. 


3, Z Procedure for Starting a Nonequilibrium Solution 


As mentioned above the nonequilibrium solution cannot be started in 
the reservoir. The solution is started by considering the nonequilibrium 
solution to be a perturbation about the equilibrium solution in the early stages 
of the expansion. Thus, as for the equilibrium flow, the nonequilibrium solu- 
tion begins with taking fixed temperature steps from the reservoir state. 

At each temperature the equilibrium composition, pressure, density, and 
velocity are computed in the same manner as described in Section Z. Z. In 
addition, the equilibrium values of » dT/d% * and dinp/d*, * which are 

needed in the perturbation calculation, are computed at each point. Equations 
(28) , (34), and (36) remain valid for infinite-rate equilibrium flow. However, 
the rate Equations (30) cannot be used to determine the remaining (s-c ) 
species gradients in equilibrium. These additional equations for the species 
gradients are obtained by differentiating the equation which expresses the 
condition of chemical equilibrium, i. e. Equation (5): 

Zi± d A . AL At 

** £7 T \ T J d* It 1 


I eCq Y 

* “ ~h 


f) 


dT 

d -jt 


(39) 



j = C 5 


The van't Hoff equation which can be written 


dJU K; 

tf- 

dT 






has been used in deriving Equation (39). Equation (39) may also be obtained 
by enforcing the conditions of chemical equilibrium ( ^ = 0 and d^ /dd - 0) 

in Equation (33), The result is then 
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<££. . 0 

X = /, ... , ^ 


[39a) 


Only ( 5 - c ) of Equations (39a) are linearly independent as may be seen by 
noting that the K if ^ may be written in terms of the ( s-c. ) Y * for the 

L ' i i 

equilibrium formation reactions (Equations (2)). 



(41) 


In other words, only the formation reactions need be considered for a flow 
in chemical equilibrium. If the fiij matrix is written for the formation 
reactions it is seen that 


/. . > >>. . f ° r j = i, ■ ■ ■ , c 

r *1 V 3 

A 7 *~\j - C for / =c + , >---* 5 


(42) 


thereby demonstrating the equivalence of (39) and (39a). 

In order to calculate the equilibrium values of , and 

** A 

— ■ f from Equations (28), (34), (36) and (39) the value of — — — must 

be known. The value of the area ratio at a temperature step is calculated from 
Equation (13) using the results of the equilibrium calculation for /o and U , 
and the equilibr ium-flow value of the critical mass flow, vn . In the nozzle 
flow program the area distribution is specified as a polynomial in # , the 
distance along the nozzle, taking the origin at the geometric throat. Having 
determined the area ratio at a given point in the calculation the specified 
formula for Afa) is inverted to obtain at the current computational point. 
Then, using this value of /£ , UZnh /d can be calculated from the specified 
description of the nozzle. 

As mentioned in the previous section the nozzle geometry may be speci- 
fied in either of two ways. When a simple shape^i. e. y a wedge, cone or 
hyperboloid geometry, is used A(#) is a quadratic expression in ^ . For 
the case of a simple geometry then, /t may be obtained from a given area 
using the relations 
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x O 


( 44 ) 


X 


A x + /a * - 4A 3 (At a) 

^ ’ 


. . B, 4-6, * > 0 |44a) 

ZB, 

The signs employed in the quadratic formula are determined by the behavior 
of A(x) in the converging and diverging sections. 


The other method for specifying the nozzle geometry is to use a series 
of polynomial fits to actual nozzle contours. In the first intervals upstream 
and downstream the polynomials must be quadratics. In quadratic regions 
Equations (44) or (44a) are used to invert the A( expression to find & . 

In the intervals employing higher-order polynomials a Newton-Raphs on iter- 
ation procedure is used to invert the A fa ) relation. The value of # at the 
previous step is used as an initial estimate in this procedure. Since quadratic 
expressions are used in the first intervals in both sectionsi initial estimates 
for /jC which lie in the correct interval are ensured. The correct polynomial 
to be used for A(x ) is found by comparing A with the values at the dividing 
points between intervals. 

Summarizing the discussion of the starting procedure to this point: 
the equilibrium flow solution is obtained as a function of temperature. Using 
the equilibrium calculation for p>u. and the equilibrium flow value of , 

A(tf) is calculated at each point; the specified relation for A in terms of 
is then inverted to find /£ . In turn, the calculated value of /£ is used 
to evaluate dJU,A/d<t at the current step. Having computed the equilibrium 
composition and flow properties and evaluated d/tiA/d# , the equilibrium 
values of dy^/et/jL > dT/d* » and dJUp/dy. maybe found from Equations (28), 
(34), (36), and (39). 


The values of the above derivatives are then used to obtain the pertur- 
bations of the composition and flow properties from their infinite- rate 
equilibrium values. This procedure is continued until the perturbations are 
large enough to start the numerical integration of the nonequilibrium solution. 
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The calculation of the perturbations and the size they must attain to start 
the numerical integration are discussed in ensuing the section. 

Once the perturbation calculation indicates that the nonequilibrium 
solution can be started, the perturbations in the species concentrations and 
the temperature are added to the equilibrium flow values. Then the static 
enthalpy and velocity are recomputed using Equations (12) and ( 1 5) ? r e spectively. 
The density is obtained from Equation (13) using the equilibrium -flow, critical 
mass flow and the area ratio computed in the starting procedure. The 
pressure is then evaluated from the state Equation (11b). Using the thus 
obtained approximations to the nonequilibrium flow composition and properties 
at a point, the slopes ■ ■ ^ , — , and J^P can be calculated from 

a. % at CL 

the nonequilibrium flow equations as discussed in Section 3. 1. Notice that 
the nonequilibrium solution then proceeds in steps of *£ instead of T as in 
the starting procedure. Also, at each new value of # the area ratio, A(&) > 
and dlnA/d# can be calculated directly from the specified relations for the 
nozzle geometry. The rest of the nonequilibrium solution then proceeds as 
described in Section 3. 1. 

If the departure of the flow from equilibrium is expected to be small 
upstream of the nozzle throat then the starting procedure may be initiated 
at the nozzle throat. Provision is made for taking a variable temperature 
step from the equilibrium throat conditions. After this initial step the 
starting procedure continues at the usual temperature interval employed in 
the frozen and equilibrium computations ( 0,01 Tj ). If the starting procedure 
is initiated at the reservoir and the solution approaches the nozzle throat 
without the perturbations becoming large, the calculation is continued in the 
same manner as if the option for starting at the throat were exercised. 

3. 3 Modification of Nonequilibrium Calculation when Integration is Started 

Upstream of Nozzle Throat 

In discussing the nonequilibrium calculation and the starting procedure, 
the critical mass flow was assumed to be the equilibrium-flow value. If the 
flow is appreciably out of equilibrium upstream of the nozzle throat then the 
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critical mass flow is not known beforehand. The results of earlier studies 
indicated that the nonequilibrium density distribution upstream of the throat 
and the critical mass flow differed only slightly from the equilibrium flow 
values. On this basis an inverse procedure is used in the upstream region. 

When the perturbation calculation (Section 3. 4) indicates that the 
numerical integration of the nonequilibrium solution should be started up- 
stream, the following procedure is adopted. The equilibrium flow density 
distribution and critical mass flow are considered to be the prescribed 
boundary condition instead of the specified nozzle geometry. The equilibrium 
density is assumed to be related to the specified area distribution by 

(/>*) (/-/**) = C ( 45 ) 

This is the form of the relation between p and A for an ideal gas where 
oc = /- / . In using Equation (45) to describe the (o(A) relation, the con- 

... it 

stants oc and £ are evaluated from the following conditions at A = 1: p “ p 
(computed for equilibrium flow) and dA/d p - 0. Differentiating (45) and 
setting A = It ciA/dp = 0, and p = p* yields 

(pf (oc +2) =2 (46) 

This resulting equation for oc is solved using a Newton-Raphson-type 
iteration to satisfy the condition 

f(oc) (pC* Z) - z = 0 (47) 


Since &(oc) must be zero, the 7i estimate for oc is given by 

„ t<*"> 


71+1 

oc ^ oc — 


3 oc 


(48) 


oC*o i 


where is obtained from a two-term Taylor series expansion about 


OC 


If Equation (46) is rewritten as 

OC Jrfv p* + JLro (c*c + 2 ) = 2 
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and £ fv, ( oc + Z) is expanded in a Taylor series about oc = 2 , retaining the 
first three terms, the following approximation is obtained for oc : 

« » 4 + 8 1™ p * < 4, > 

This result is used as an initial estimate for oc in the iterative solution of 
Equation (46). Once oc is evaluated then C may be found by substituting 
p = p and \ in Equation (45). 

Having found the constants o c and C > the density-area relation may 
then be employed in the nonequilibrium solution upstream of the throat. At 
each new value of 4 , in the numerical integration the area ratio is found from 
the specified geometry. The equilibrium flow density and are then 

evaluated as follows. Again a Newton-Raphs on- type iteration procedure is 
employed to find the value of p which satisfies Equation (45) for the given 
values of A , oc , and C . In this case the relation to be satisfied is 




(50) 


and the 71 + / value for p is obtained from a three-term Taylor series 
expansion about the n ^ estimate 


* r(r> 
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= 0 


This equation is then a quadratic equation for */° = /° ~ p 

is given by 
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D, 


(51) 


The solution 


(52) 


where 

v--f(p v ) 


J>,« 


df 


P*P 


-n A£ 

A = zp* 


~ ft 


P--P 


The sign in the quadratic formula is chosen so that the smaller M is used. 
The density from the previous computational step is used as an initial estimate 
in solving for p at the new . 
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In order to compute the nonequilibrium solution when p (j)C ) is 
specified dirt p/dtf must be calculated. Once the equilibrium flow density 
has been evaluated, dln,p/di£ may be calculated using the values of A and 
dJLb/cU corresponding to the actual specified geometry, and the relation 

_ 2 c aU^A 


d# 


o <( / oA)*-(oc + 2)C d# 


(53) 


This relation may be derived by logarithmically differentiating Equation (45). 
The nonequilibrium flow Equations (28), (30), (34) and (36) are then solved 
for df^d#* dT^d# » and dl^A/d'jC > where A is the cross-section of the 
nonequilibrium streamtube flow having the density variation pOt) , and the 
numerical solution proceeds. For a converging-diverging nozzle d ZrvA/d% 
passes through zero at the geometric throat. In the governing equations for 
the flow of a reacting gas through a nozzle, a branch point exists where the 
flow velocity becomes equal to the ’’frozen” sound speed. Since the frozen 
sound speed for a reacting gas is always slightly larger than the velocity at 
the point of minimum cross section [i. this branch point will 

occur slightly downstream of the geometric throat. ^ Hence the solution for 
the upstream region is continued until the value of the area ratio, A (d) , 
reaches a value which is typically downstream of this singular point in the 
solution. Then the solution is continued in the usual way, i, e. by specifying 
A fa) and solving for p ( In the downstream region the specified A 

-V 

is modified to obtain a smooth transition from the am obtained using the 
inverse procedure in the upstream region. 

When the nonequilibrium calculation is switched from the upstream 
region to the downstream region the calculated value of d/tvA^/dt is com- 
pared with the value given by the modified area-ratio relation for the down- 
stream section. If these two values differ by greater than a specified amount, 
the switch from the upstream to downstream region is made at a smaller 
value of d . Also, once the switch from the upstream region to the down- 
stream region has been effected ? the value of dtr^p/d'A computed at each point 
in the nonequilibrium solution is monitored. If dim*p/diC becomes positive 
slightly downstream of the switching point from the upstream to the downstream 
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region, this switching point is moved further downstream. If cL£*,p/dt be- 
comes positive at a significant distance from the upstream to downstream 
switching point the reason cannot be attributed to the branch point in the 
solution and the calculation is stopped. 

As evident in the present discussion, when the solution is started in 
the upstream region the area-ratio variation of the nonequilibrium solution 
is slightly different from the specified geometry. Consequently, the area- 
ratio variation for this case does not become unity at the throat. In principle 
this area ratio distribution should then be renormalized using the value at 
the minimum cross-section. In practice, this value differs from unity by 
less than 1% and renormalization is unnecessary. 


3. 4 Perturbation Calculation Used in Starting Nonequilibrium Flow Solution 

As mentioned earlier the numerical integration cannot be started at 
the reservoir because there the gradients of the flow properties vanish. The 
procedure adopted in the present program is to start the numerical integration 
by considering the nonequilibrium flow solution to be a perturbation about the 
infinite -rate equilibrium flow in the initial stage of the expansion. Let the 
species concentrations and flow properties be given by 

y. = +8r- 

T = T + ST (54) 

p = p + Sp 

where the bar symbol denotes the infinite rate equilibrium value. The element 
conservation equation (Eq. (27) ) yields 

5 

°V* ^ - 0 -*. = », 2, c (55) 

Similarly substituting (54) in the species conservation equations (Eq. (3 Oy ) 
yields, to first order, 


d *i , 

eC % d-v 



Pxl l + ’l,S p ;'-P; 
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The bar 


Now an important point in the perturbation procedure arises, 
symbol denotes that a quantity is to be evaluated using the equilibrium flow 
values of the dependent variables. Thus, R does not denote the value of 
for the infinite - rate equilibrium flow, but rather Ft evaluated using the actual 
value of Jt and p , f , and /■ . Since, £ . is the same as in the equilibrium 
flow, = 0 and the above equation becomes 



The 8l’* may be written in terms of S /■ > 8T > and retaining only first- 
** Z 3 ' 

order terms as 


Sr st + s r +X 

ar dp r yrr 




Evaluating the coefficients of the perturbations leads to 

<5 




1,9 


I, = Sy rt 4 £ (¥- i ) st '£- S/ ’ 


/•» 7 

From the energy Equation (15) 

3 H + U 8 u - 0 

Using Equations (12) and (13) this becomes 


%t* { ST -*£- 

U 1 f: - 1 3 * 77i 2 dT F 




- o 


(57) 


(58) 


The remaining governing equation, the momentum equation, has not as 
yet been perturbed. The perturbed momentum equation would involve 
derivatives of the perturbations. The use of the perturbed momentum 
equation would then make the solution for the perturbation quantities quite 
difficult. At this point two approximations based on the pure diatomic gas 
studies^ are introduced 


and 


8 p - 0 


(59) 

(60) 
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In Reference 9 it is shown that the series expansion procedure employed to 
start the solutions obtained in Ref. 7 justifies these assumptions near the 
reservoir. 

Also discussed in Reference 9 are alternative methods which might 
be used instead of Equations (59) and (60). Rather than using 6^ = 0 in place 
of the perturbed momentum equation the condition 

8S' r 0 (61) 


was employed. Based on numerical computations with the program Equation 
(61) is known to be accurate to higher order than &p - 0* Equation (61) is 
used in the program modification discussed in Appendix A, Also, a successive 
approximation scheme was employed to estimate the terms d($ d in 
Equation (56). This improvement offered no significant advantage over the 
method presented here. 


Using Equations (57) and (60) in (56) yields 



ft s 







(62) 


Now Equations (55), (58), (59), and (62) provide ( s + 2) algebraic equations 
for <5/. i ST > and 8p in terms of the equilibrium flow properties at a given 
nozzle station. 


As described in Section 3. 2 the nonequilibrium nozzle flow solution is 
started by obtaining the equilibrium flow solution at successive steps in tem- 
perature. At each step the results of the equilibrium calculation are used 
to obtain the values of <5 ?.■ , St , and So . Then, using these results, values 
of 6 1 . are computed at each point from Equation (57). When any one of the 
exceeds a specified value, the values of <5/. , §T > and Sf> are added 
to the respective equilibrium flow values and the numerical integration is 
begun. 


Numerical difficulties arise when trying to integrate the nonequilibrium 
equations for small values of ^ . , Consequently, the solution should be 

started for as large values of 6* af , i. e., as far downstream as possible 
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without the linear perturbation procedure becoming inaccurate. At each 
computational step in the starting procedure the perturbations are added to 
the equilibrium values of , T, and and values of £ , are computed 

from Equation (33). Since ;jr = 0, if the linearization procedure is valid 


(?• + 8 /■ , T * ST, P + 8p) & 5 1 

0 / 


(63) 


Experience obtained in the use of the program indicates that Equation (63) 
holds for values of 6 £ . up to the order of 10 Thus, the numerical 
integration is started when the largest § 7 is typically 0. 1. 

In the above discussion it is indicated that all of the species gradients 
are evaluated from the rate equations once any 8i. exceeds the specified 
test. In the nonequilibrium flow of a multicomponent mixture it is possible 
for some reactions to remain much closer to equilibrium than others. In 
Reference 9, a procedure is described whereby the numerical integration is 
started at different points for different species. However, if there are more 
reactions than species there is more than one production term in the species 
conservation equations and the procedure becomes inconvenient. Further- 
more, the method does not help when two competing reactions keep a species 
concentration near the equilibrium value. As discussed in the ensuing section 
the numerical-integration scheme employed in the computer program was 
developed specifically for strongly coupled nonequilibrium processes. The 
use of this method has alleviated difficulties associated with some species 
remaining near their equilibrium distributions. 


3, 5 Numerical Integration Method Employed in Nozzle Flow Program 


The solution of the nonequilibrium nozzle flow problem requires the 
numerical integration of a system of first-order differential equations of 
the form 




(64) 


In the original version of the program a fourth-order Runge-Kutta scheme was 
employed. For this method, the formula for each in an interval is of 

form 
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where 


A A * 

' 6 

+ *f 4 ) 

(65) 

, A * 

*•=*'* 2 

y. ■ 7' * *. V 


A 

*' + z 

fs ’ 7, * 

(66) 

= /£, 4- 




and 

£• = f fe > ^ ) 

The 4th order Runge-Kutta formula (Equation (65)) is derived assuming that 
the slope in any integration interval is given by 

A + + (*-*,)* (67) 

The Runge-Kutta method is generally quite adequate for the numerical 
solution of nonequilibrium flows. However, when the flow is near an 
equilibrium condition this method becomes unstable except for extremely 
small step sizes. The way in which the instability arises can be explained 
in the following manner. Near an equilibrium condition the species conserva- 
tion equations can be written in the form 

d< 


d 4, 


t- S' -p(y-y) 


( 68 ) 


where p , an inverse relaxation length, is a large number and y^, 

the departure from equilibrium, is small. Differentiation of Equation (68) 
shows successive derivatives of y to be p times larger than the previous 
order, hence the higher order derivatives of ^ are large. 

The 4th order Runge-Kutta method is inadequate for such equations, 
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Recently, Treanor has developed a modification of the 4th order Runge- 
Kutta scheme which alleviates the numerical inaccuracy of the original 
method for near -equilibrium flows. Instead of Equation (67) the following 
relation is used for the slope in an interval. 
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(69) 
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This relation leads to the following integration formula 


Ay = Ad |a F ( t B (A d ) F t + C (Ad) F 3 'J 


(7 0) 


where 
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A= f, 

8= A (72 , 

c * (kt ) 1 ^ * 1 p y.) ' fc ^ 1 p y.) ' - fc * p *. HV p y, )] 

The f. and are again given by Equation (66). The modified Runge- 

Kutta expression is of higher order than the 4th order Runge-Kutta expres- 
sion but approaches the latter for small values of PM- 


When P( Atf)is large the value obtained for y + using the Runge-Kutta 

relation (Equation (66)) is inaccurate. Since P is evaluated knowing f ( , 

■P and p then u can be evaluated using Equation (67) with the quadratic 
T » * 

term omitted. The result is 

% = y l +&*fef,F.+f,fc-2Fj + f»F x Pd-*)} (73) 

When is small, which corresponds to a situation where 

Equation (68) is not valid, negative values of P are possible. In this case 
the original Runge-Kutta formula is used. Furthermore, when , 

and f are nearly equal } P cannot be evaluated accurately using Equation (7 2), 
However, when the slopes ^ , and are very close, the Runge-Kutta 

formula should be accurate and again it is used to obtain Ay . 
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In a nozzle flow calculation the numerical integration is stable for a 
much larger step size in the later stage of the expansion than in the initial 
stage. The numerical integration is begun with a step size which is known 
to be small enough to ensure stability. Then if the solution passes certain 
integration tests the step size is increased. In this way much larger inte- 
gration steps are taken in the downstream part of the solution and the com- 
puting time is not prohibitive. 

Two integration tests are used to control the step size in the program: 

(1) the percentage changes of and T in an interval must be less than a 

$ 

specified value (2) and T must be positive at both intermediate points 

and the final point in the integration. Initially, if these tests are passed 
on a specified number of successful steps then the interval size is increased 
by a factor of 1. 1. If the interval size has been increased successfully a 
specified number of times then the factor by which the interval size is 
changed is increased by 0. 1. In the program the number of successful 
increases required before increasing the rate of increase is taken to be the 
same as the number of successful steps required before increasing the 
interval. When any of the integration tests are failed the step size and the 
rate of increase are reduced. 

The use of the modified Runge-Kutta scheme including the latter 
method for determining ^Equation (7 3)^ has considerably extended the 

capability of the nozzle flow computer program. Many solutions have been 
obtained which could not be obtained using the 4th order Runge-Kutta 
method. Also, solutions can be obtained with far fewer integration steps 
and hence for smaller computing times. The results obtained for the species 
concentrations in a solution for a nozzle expansion of air from reservoir 
conditions of T^ = 10,000 °K ; 'p/ c = 4400 atm. are shown in Figure 1. As 
discussed in Reference 11 these reservoir conditions fall in the range where 
the nitric-oxide shuffle reactions keep the nitrogen-atom concentration very 
near its equilibrium value throughout the expansion. This solution, which 
was originally obtained using the Runge-Kutta method, was repeated using 
the modified Runge-Kutta scheme. A comparison between the two results 
for the nitrogen-atom concentration for a section in the initial part of the 
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expansion is shown in Figure 2. As can be seen the step size in the modified 
Runge-Kutta method is larger. Attempts at increasing the step size with 
the Runge-Kutta method lead to oscillations of the character demonstrated 
in Figure 2. 

In the above comparison of the modified and 4th order Runge-Kutta 
schemes the Integration tests which limit the relative changes in the species 
concentrations and temperature were set at 1%, The modified Runge-Kutta 
method can be used successfully with larger integration tests: the integration 
tests on the species concentrations and the temperature arc typically set 
at , 10 and , 05 respectively. When these values of the integration tests arc 
used, solutions can be started with integration steps as high as 10 cm for 
the modified Runge-Kutta scheme as opposed to 10 ^ cm for the 4th order 
method. 
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4. THERMODYNAMIC AND CHEMIC AL-KINETIC MODELS 


In this section the methods employed in the computer program for 
specifying the thermodynamic and chemical character of a gas model are 
described. As mentioned in Section 2. 1 the chemical species which are to 
be considered in a mixture are specified through the . matrix. Further 
data which must be supplied to compute the thermodynamic properties of the 
mixture at a given point are thermodynamic properties » zt* , 

o / 

and 5* for each species. In the nonequilibrium computation, data for 
chemical reactions among the species must be provided in addition to the 
thermodynamic properties of the species. 


As mentioned previously, it is assumed that the expanding gas mix- 
ture remains in thermal equilibrium and consequently that there is no 
coupling between the chemical kinetics and the internal degrees of freedom 
of the gas. In high temperature expansions vibration-dissociation coupling 
may become significant. Unfortunately there are not at present any estab- 
lished models available for the assessment of these effects, hence they have 
not been included in the program. 

The data for specific kinetic models are not included in the present 
discussion. However, the reports and papers describing the applications of 
the program* 1 * contain detailed summaries of the models employed. A 

sample kinetic model for air and the preparation of input data for this model 
are illustrated in Section 6. 2, 


4. 1 Specification of Thermodynamic Properties of the Species 


For a mixture of ideal gases the thermodynamic properties , 

o f 

and s . are only functions of temperature. Two methods are employed in the 

o 

computer program to specify these required functions of temperature (a) a 
simple - harmonic -oscillator model, (b) polynomial fits to more accurate 
calculations of species properties. In both methods the species are assumed 
to be vibrationally and electronically excited, in equilibrium with the local 
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translational temperature. Both the translational and rotational degrees of 
freedom are assumed to always be fully excited. In the harmonic -oscillator 
description the electronic excitation is assumed to be given by the first few 
terms in the partition function sum. In the thermo-fit method, as the poly- 
nomial fit method is called, the computations of the vibrational and electronic 
contributions are typically more accurate than the simple harmonic oscillator 
calculation since account is usually made of anharmonic vibrations, vibration- 
rotation coupling, and more terms in the electronic partition function. 


Using the simple harmonic -oscillator model the specific enthalpy for 
a monatomic or diatomic species may be written 

-v r 


A io = 5 + 2fo,-i; 
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i=i i 


( 74 ) 


where >1 . = no. of atoms per molecule of V th species, 

J 

Q = characteristic vibrational temperature of J th species, 

i / 

= degeneracy and energy of i th electronic state of j, th species, 

l i i 

and . is the formation enthalpy of each species at standard conditions 
(usually zero temperature). The maximum value of Lj , the number of 
terms in the electronic contribution to A; . allowed for in the program 
is 10. 


For the oscillator model the chemical potential at standard pressure 
(1 atm. ) is given by 
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and 


A, = 4 U ~(n r \)l*Q' 

p' / </ 


(77) 


'3 Z A 

In Equation (7 5) 7nj , P , ,^j/'are respectively the mass of the ^ th 
particle, Boltzmann’s constant, Planck's constant, and the standard 
pressure, all in c.g.s. units. Also is the characteristic rotational 

temperature in °K, 

In the nozzle flow computations the quantities Sj and are also 

needed. Using Equation (17) 5^ may be found once and yt-y are known, 

(17) 


7 T 


Differentiating Equation (74) 
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Thus, in the harmonic -oscillator description of the species properties the 


following constants must be specified for each species: 
0 


V 


and 


n i * 


J 




L i ’ \ 

The thermo-fit method of specifying or yj_ . is to fit a fourth- 

/ a 

order polynomial in temperature to calculated data for the individual 


species. Thus 


o v 

* T - = *i -^ T ')-^ f T - %(T‘) - (79) 


A, is obtained by differentiation 

t 




= a.j + P- t'+ c.(j'S +dj(r') + e-(j'f 


(80) 
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Again 5* is computed from Equation (17) and C is obtained by 

i . * 

differentiating Equation (80). Thus, in order to use the thermo-fit method 


the polynomial coefficients cl* , Jr; , c- , ci / 

0 4 i oo 

with Ji . must be specified for each species. 


and together 


In the computer program provision is made for using the harmonic- 
oscillator description for some species and the thermo-fit method for others, 
in a given calculation. The use of this option is discussed in Section 6. This 
option is useful, for example, for a mixture containing triatomic or more 
complex species since the above harmonic -oscillator description can only 
be used for monatomic and diatomic species. 


Another option incorporated in the program allows the use of a thermo- 
fit description over the initial part of the expansion and a harmonic -oscillator 
description downstream. This additional feature is useful for high-expansion 
nozzle flows at high reservoir temperatures. In general the harmonic- 
oscillator model becomes inaccurate at high temperatures. Thus a 
polynomial fit to more accurate data can be used in the initial, high-temper- 
ature portion of the expansion. However, the temperature variation in such 
an expansion is large and hence may decrease below the lower limit of the 
range over which the polynomial coefficients are valid. In this case, the 
harmonic - oscillator description which is accurate at moderate and low 
temperatures (up to about 5000°K in air) can be used to extend the thermo- 
dynamic description to cover the entire expansion. The use of this latter 
option is illustrated for the air flow example presented in Section 6. 


4. 2 Chemical Kinetic Model 


The chemical reactions among the species in the gas model are 
specified through the P* • and ■ matrices. As seen from Equation (26) 

0 9 

these matrix elements represent the stoichiometric coefficients of the j th 
species on the reactants and products side, respectively, of the j. th reaction. 
In addition, the reaction rate constants must be specified for the forward 
direction of each reaction. Thus each reaction should be written so that the 
rate constant is known for the forward direction. The reverse rate constant 
is then computed through the equilibrium constant. 
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In the computer program forward rate constants are written in the form 


4 
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(81) 


which is the usual form in which experimental data are given. If third-bodies 
appear in dissociation-recombination reactions and the rate -constants are 
the same for some of the third bodies the following simplification is possible. 
For reactions having the same rate constant but different third bodies only 
one reaction is written, omitting the third bodies in the j) • * and 9. . matrices. 
Then the forward rate constant for this reaction is written 


■Af, = A :T 






^ y u - *<- 

A t > \ 


(82) 


where U^. = 1 for all species which are third bodies in the Ji th reaction and 
zero otherwise, As can be seen from Equations (30)-(33) when = 1, this 
procedure for reducing the number of reactions in the kinetic model is simply 
a regrouping of the chemical production terms in the species conservation 
equations. If there are no third bodies in a given reaction then = 0 and 
no entry is necessary in the corresponding row of the (J . . matrix. 

When the numerical solution is begun at a point where the flow is out 
of equilibrium and hence all the flow properties and the composition are 
specified^ there is no requirement for a minimum number of reactions in the 
kinetic model. However, when the perturbation method is employed to 
start from an equilibrium state, there must be as many linearly independent 
reactions in the chemical kinetic model as there are formation reactions in 
the equilibrium calculation. In other words, a sufficient number of linearly 
independent reactions must be included so that the rank of the ^ , matrix 
is (S-C) . If the rank of A . . is less than (s-c ), the perturbations in species 
concentrations are not uniquely determined by the perturbed rate equations. 

In each nozzle flow calculation the rank of ^ t , is tested using a standard 
diagonalization procedure. If the reactions thought to be important in a 
calculation do not have a /$ . . of rank (s-c), appropriate formation reactions 
with fictitiously low rate constants should be included in the model. 
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5. ADDITIONAL CAPABILITIES OF NOZZLE FLOW PROGRAM 


The nozzle flow computer program contains options which permit 
including ionized species and nonequilibrium ionization reactions, vibrational 
excitation frozen at the reservoir value, and effects of moderate gas 
imperfections on the initial, near - equilibrium portion of the expansion. The 
analysis underlying these options and the procedure for exercising them are 
discussed in this section. 

5. 1 Nozzle Flow Calculation for Ionized Flows 


The procedure for computing the flow of a reacting mixture through a 

converging -diverging nozzle can be extended to ionized mixtures by treating 

15 

electrons as a separate chemical element . Each of the positively 
ionized species would then be formed by subtracting electrons from the 
corresponding neutral species; negatively, by adding. The conservation of 
elements equation (Eq. (27)) for the electrons then enforces conservation of 
charge. The formation enthalpy of each of the ions is the sum of the 
formation enthalpy and the ionization potential of the corresponding neutral 
species. Accordingly the electrons are referenced to the same zero level 
of energy as the molecular species. The sample calculation described in 
Section 6 includes ionized species. 


In the reservoir and equilibrium-flow computations the electrons are an 

independent species or component. Hence, the concentrations of the ions in 

the mixture are found in terms of the electron concentration. In an equilibrium 

expansion to low temperature the electron concentration and thus the ion 

concentrations approach zero. Thus, when the electron concentration 

decreases below a specified value, the ionized species are dropped from the 

computation. In order to perform this truncation of the model the electrons 

and ions must be listed in a certain order in the oc L . matrix. The 

$ 

instructions for ionized flows are included in the description of the input-data 

format in Section 6. At present the ionized species are dropped from the 

-30 

computation when the electron concentration becomes less than 10 moles 
per gram of mixture. However, the electron concentration cannot be 
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expected to be calculated accurately at these very low levels. The accuracy 
of the equilibrium solution for an ionized flow is further discussed in Section 
7 under the description of the subroutine used for solving the set of linear 
algebraic equations arising in the Newton-Raphson iteration. 

If the electrons are treated as a separate chemical species the pro- 
cedure for specifying reactions between ions and neutrals or electrons, ions, 
and neutrals is the same as for chemical reactions of neutral species. 

Solutions for nozzle flows of high-temperature air have been obtained 

1 5 

including the effects of nonequilibrium ionization . The estimates for the 

rate constants of three-body deionization paths included in the computational 

model led to values of A ^ (Eq. (82)) which exceed the digital capacity of the 
3 8 

machine ( > 10 ). Since the product appears as a multiplying 

•A, 

factor in each term of the species production terms (Eq. (30)) the following 

remedy was adopted. The constant factors in all of the rate constants were 

3 8 

lowered by a factor such that the largest value would be -^10 , Then the 

characteristic length, JL , was increased by this factor. When this 
procedure is employed the coefficients in the Afa) relation must be 
compensated so that the correct value of A ‘/A' is given at each s 

In previous experience with this option of the program it was found that 
when several positive ions were included in the calculation, charge was not 
conserved in the nonequilibrium solution. This discrepancy has been 
attributed to integration error arising from using the differentiated form of the 
element conservation equations (Eq. (28) ). This problem does not arise in 
conserving the neutral elements, element conservation being obtained to 
within four-decimal places in existing solutions. The lack of charge 
conservation has been remedied by equating the electron concentration to the 
algebraic sum of the ion concentrations after each numerical integration step 
in the nonequilibrium solution. 

Using the above procedure for solving ionized flows is valid only for 
weakly ionized flows in the absence of applied fields. At very high temper- 
ature and low pressures gases become highly ionized and the assumption of 
a mixture of perfect gases breaks down. Attractive forces between the 
charged particles become important and Debye shielding effects must be 
considered in the state equation. Also the present analysis neglects 
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diffusion and nonequilibrium electron temperatures, assumptions which are 
not valid for highly ionized nozzle flows. 

5. 2 Nozzle Flow Calculation for Frozen Vibrational Excitation 


In describing the nozzle-flow computation for frozen, equilibrium and 
nonequilibrium chemistry the vibrational degree of freedom was said to 
remain in equilibrium with translation during the expansion. The other 
limiting assumption for the degree of vibrational excitation, i.e. frozen at 
the reservoir value, is also allowed for in the nozzle-flow program. This 
option may only be exercised when the harmonic -oscillator model is used 
for the thermodynamic data because that model explicitly displays the 
vibrational contribution to the internal energy. 


If the vibrational excitation is assumed to remain frozen at the 
reservoir value, the equation for the specific enthalpy of the ^ th species 
becomes 
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Accordingly, the chemical potential of each species is given by 

-{v - < 75 
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The derivative now contains no vibrational contribution and hence is 
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given by Equation (78) with the vibrational term omitted. s . 
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5. 3 Correction for Moderate Gas Imperfections 


At very high densities, the assumption that the molecules of a gas are 
noninteracting, point particles is no longer valid. At high temperatures and 
high densities, if ionization is negligible, the dominant intermolecular 
forces in a gas mixture are the repulsive forces between particles. Thus, 
in this situation the gas imperfections essentially stem from the finite volume 
of the particles . 


When these effects are moderate, i.e. when the volume of the particles 

is 10% or less of the volume occupied by the gas, they can be approximated 

by assuming the mixture to be composed of fictitious hard-sphere particles. 

2 1 

The equation of state for such a mixture is 


r = 


p 


i - 


££ o 



(11c) 


where is four times the volume of the particles, is the molecular 

weight of the undissociated mixture, and ^ is the molecular weight of the 
mixture. The so-called "co-volume" of the particles, Jr 0 , is given by 


J- 

O 


IT N a cr 


(83) 


where N 0 is Avogadro's number and cr is the diameter of the particles. 
The state equation may be written 


■p' zj'/n 


(84) 


where jd ' is an effective density defined by Equations (11c) and (84). 


When a hard-sphere model is employed to account for gas imperfections 
there is no correction to the ideal-gas internal energy. However the enthalpy 
is given by 


h'= e'+ p/p' = 


?' *.T' 
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or 


H 
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(85a) 


If the temperature and pressure of a gas mixture are specified and the 
diameter of a hard sphere representing the mixture is used, then the gas 
density is given by (11c) and the enthalpy is given by (85). 


Using the above approximation the effect of gas imperfections on the 
nonequilibrium nozzle flow calculations can be accounted for in the following 
manner. The reservoir computation is carried out for a specified 
temperature and pressure. The composition is the same as for an ideal 
gas mixture but the enthalpy and density are computed using (85) and (11c). 
The computation for the equilibrium portion of the expansion proceeds in the 
same fashion as before. However, now the thermal and caloric state 
equations are written 

i ^ \ 
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The solutions available for nonequilibrium nozzle flows indicate that 
high reservoir densities delay the departure from equilibrium. At densities 
sufficiently high for gas imperfections to be significant, the flow will remain 
near equilibrium over an appreciable portion of the expansion. Thus, in the 
nozzle-flow program the corrections for gas imperfections are included in 
the state equations only for the near - equilibrium portion of the flow. The 
flow will most likely have expanded to the point where the ideal-gas 
equations apply before departing from equilibrium. 
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The above procedure for accounting for gas imperfections has been 

applied to nozzle flows of air for reservoir temperatures of 8000, 10, 000, 

12, 000 and 15, 000°K and reservoir densities of 10 and 100 times sea level 
1122 

density. Hilsenrath has computed the thermodynamic properties and 
composition of equilibrium air including second virial corrections for 
temperatures up to 15000°K densities up to log ^ /°fp = 2.2. Considering 
air to be composed of fictitious n air molecules”, a value of or of 2.6Ain 
Equation (11c) was found to fit Hilsenrath’s calculations for the compressibility 
factor 2 very well in the range of interest. For given temperature and 
pressure, the correction to the density at an ideal gas value of 100 f>° is 
about 10% in this temperature range. Also, at this density level Hilsenrath's 
results indicate the corrections to the internal energy and enthalpy to be 
about 1% and 3% respectively. While the model used in the nozzle flow 
program assumes no correction to the internal energy, the correction to 
the enthalpy given by Equation (85) is about 3%. In the nozzle flow solutions 
obtained in Reference 11, the gas had expanded to the point where it could 
be considered ideal before significantly departing from equilibrium. Further 
details on the application of this method to high density airflows and the 
results of such calculations are presented in Reference 11. 

Thus, the computer program is capable of including the effects of 
moderate gas imperfections on nozzle expansions for high reservoir densities. 
In order to use this option an empirical value of O' , the diameter of a 
"molecule” of the mixture, must be chosen for the conditions of interest. 

As discussed in Section 6, this option can be excluded from the calculation 
by simply setting J- =0. Then the equations in the program reduce to 
those for a mixture of ideal gases. 
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6. INPUT-OUTPUT PROCEDURES 


This section describes the input-output procedures for the nozzle-flow 
program. The preparation of input cards is discussed, including a description 
of various options available in the program. This is followed by sample input 
data for an ionized air calculation. The output format is then described, 
illustrated by sample output data for the ionized air example. 


6. 1 Preparation of Input Cards 

The input cards necessary to initiate a nozzle flow computation are 
listed below. The cards are listed in the order they appear in the data deck 
and an explanation of the information required is given for each card type. 
The definitions of the FORTRAN symbols used to indicate the information on 
each card type are given in Section 7*3; sample input cards are presented in 
Section 6. 2. 


Card T ype 1 
ISW1A 


ISW2A 


ISW3A 


ISW4A 


ISW5A 


Format ( 1614) 

1 if solution with frozen chemistry is required 

0 otherwise 

1 if solution with nonequilibrium chemistry is 
required 

0 otherwise 

1 if solution with equilibrium chemistry is 
required 

0 otherwise 

1 if another run follows 
0 otherwise 

0 this indicator is not currently used 
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ISW6A 

ISW1B 
IS Vv r 2 B 
ISW3B 
ISW4B 
ISW5B 
ISW6B 

Card T ype 2 
ACOM(I) 

Card Type 3 

IRUN 

ISC 

ISS 

ISR 

IC 

NQS 

IUPD 


1 if reservoir calculation only is required 
0 otherwise 

Six extra indicators are available for additional 
options as required. Storage is allcted in COMMON 
for those indicators. If not used, the remainder of 
the card may be left blank. 


Format ( 12A6) 


up to 120 alpha-numeric characters describing 
run, printed as heading on first page of output, 
(There will be 2 cards.) 

Format ( 1614) 


run number 

number of chemical elements (maximum of 10, 
including the electron) 

number of chemical species (maximum of 20) 

number of reactions (maximum of 64) 

number of ions (do not count the electron) 

number of successful integration steps before 
increasing step size (normally 4) 

1 if starting upstream of throat 

0 if starting downstream of throat 

The nonequilibrium solution may be started either 
upstream or downstream of the nozzle throat. 

For upstream starts, the solution is begun by 
stepping off in temperature an amount DELT3 
from the reservoir. For downstream starts the 
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NFIT 


KHO 


KKUR 


NAFIT 


INEQV 

IROBAR 


Card Type 4 

CTAP 

PRESA 

CTMXX 


program computes the equilibrium conditions 
at the throat and then steps off in temperature 
DELT2 to begin the nonequilibrium solution. 

The downstream start may be used for cases in 
which nonequilibrium effects are negligible 
upstream. However, if any question exists the 
solution should be started upstream. 

1 if any thermo-fit cards are to be used 

0 otherwise 

1 if any harmonic -oscillator cards are to be used 

0 otherwise 

1 if third body matrix is used (see Card Type 8, 

QQ(D ) 

0 otherwise 

1 if fitted nozzle geometry is used 

0 if standard nozzle geometry is used 
(see Card Type 13a, b) 

1 if frozen vibration 

0 if equilibrium vibration 

1 if effective density, RHOBAR, is printed 
(see Section 6.3) 

0 otherwise 

Format (7F 10.0) 

reservoir temperature in °K 
reservoir pressure in atm 

temperature at which switch is made from 
thermo-fit to harmonic-oscillator data (in °K) 

(see Card Types 14-18) 
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CXMAX 


SL 


B ZERO 


T ST OP 


Maximum value of <£ 

The nonequilibrium calculation may be terminated 
either by an upper limit on ^ or a lower limit on 
temperature (TSTOP). Both values must be read 
in, the one not being used set to a large (small) 
number. The frozen and equilibrium solutions 
use only the temperature stop. For these 
calculations TSTOP must be greater than the size 
of the temperature step used (DELT1, which is 
set equal to .01 T/ in the program) otherwise a 
calculation will be attempted at nearly zero 
temperature . 

Characteristic length used for nondimensionalizing 
(cm.). When a standard geometry is used the 
nozzle scale may be changed by simply changing 
l . The equations for a wedge, cone, or 
hyperbolic nozzle become A = 1 + ^ , 

A = 1 + 2 -f , and A = 1 + < if jt is 0 , 

and A,* /ta*s V , respectively. In these 
expressions for X , X is the half-height of the 
wedge throat, a * the throat radius of the conical 
or hyperbolic nozzle, B is the half-angle of the 
wedge or cone, and is the half-angle of the 
asymptote cone for the hyperbolic, axisymmetric 
nozzle. For fitted geometries Ji is usually set 
at 1 cm so that the coefficients in the area-ratio 
expression are the same as for ^ measured in 
cm. 

Constant ty* used in imperfect gas correction. 
(Section 5.3) (cm^/mole) set at zero if ideal gas 
equation of state is used. 

temperature stop in °K 
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Card Type 5 
T PRINT 

DELTAX 

DELT2 


DELT3 

TTEST 

GTEST 

Card T ype 6 
CECHII(I) 


Card Type 7 
ELMENT (K) 


Format (7F I 0. 0) 

Temperature interval between printouts 
(nondimensional) if TPRINT = 0, every step is 
printed . 

Initial integration step size (nondimensional), 
usually 0.01. 

Temperature step from throat for starting the 
nonequilibrium solution downstream (non- 
dimensional), usually 0.01. DELT2 must be 
small enough so that departures from equilibrium 
are still small at the starting point. 

Temperature step from reservoir for upstream 
start, usually .001. 

Integration test: maximum size of 
Runge-Kutta step, usually .05. 

Integration test: maximum size of 
one step, usually 0. 1. 

Format (7F 10. 0) 

The switch from equilibrium to nonequilibrium is 

made when 8 lies in the range CECHII(I) 

£ S l. ^ PCTEST (CECHII(I)). If S / 

crosses this interval in one step, the step size 

is reduced until a value of 8 X „ is found in the 

t L 

correct range. The switch is made when any one 
of the S % />d/ reaches the test value. Usually 
CECHII(I) =0.1, PCTEST = 1.2. There will be 
ISR of the CECHII(I), 7 to a card. Use as many 
cards as necessary. 

Format (A6, 2F10. 0) 

Chemical symbol for element 
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CAPQ(K) 

CMW(K) 

Card Type 8 
CAI(I) 

ETAI(I) 

CEACT(I) 

QQ(I) 


Card T ype 9 
ALPIJ(I, J) 


Atoms of element Ji per molecule of mixture 
(atom fractions of elements). 

Molecular weight of element 
(There will be ISC of these cards.) 

Format (4E 14.8) 

A^, cm /mole sec(°K) * or cm°/mole sec ( # Kp 

, power of temperature dependence in Jk.^ % 

E , cal/mole °K 
act. 

A. 

1 if reaction involves a third body 
0 otherwise 

For reactions involving a nonreacting third body, 

A + B + M AB + M 

it is convenient to be able to treat all third bodies 
having the same reaction-rate constant as a single 
molecule, rather than including a separate 
reaction for each. This is made possible by 
setting QQ(I) = 1. The species which are to 
serve as third bodies are specified on Card 
Type 12. (Section 4. 2). 

(There will be ISR of these cards.) 

Format (24F3.0) 


Number of atoms of element J in species I, 

J = 1, ISC. The independent species or 

components must be listed in the first ISC rows. 

Any set of species may be chosen as components 

as long as the resulting oC^ • matrix is of rank 

o 

ISC. However, a species whose concentration 
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Card Type 10 
XNUIJP(I, J) 


Card Type 1 1 
XNUI J (I, J) 


Card Type 12 
KUR(I, J) 


Card Type 13 a 
ASUB(I) 
BSUB(I) 


may go to zero during the expansion should not 

be chosen. Thus, it is sometimes better to 

choose diatomic molecules (N^, O^) than atoms. 

If ionization is considered, the electron should 

be chosen as one of the components. Provision 

has been made to truncate the oc.. matrix if the 

^ -3 0 

concentration of electrons becomes 10 

Consequently, the electrons should be listed in 

the first row and all ions in the last rows of 

oC * * (see Section 5. 1). 

3 

(There will be ISS of these cards. ) 


Format (24F3.0) 

Number of molecules of species J on product 
side of reaction I, J = 1,. . . , ISS. 

(There will be ISR of these cards.) 


Format (24F3.0) 

Number of molecules of species J on reactant 
side of reaction I. J = 1, . . .,ISS. 

(There will be ISR of these cards.) 

Format (301 1) 

1 if species J serves as a third body in 
reaction I 

0 otherwise 

(There will be one card for each reaction which 
involves a third body.) 


Format (5E 14. 8) 

upstream area coefficients, 1 = 1, . . . . , 10 
downstream area coefficients. I = 1, ...» 31 
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5 


ATP(I) 


Card T ype 1 3 b 

ASUB(I) 

BSUB(I) 


area transfer points. 1 = 1,..., 

(These cards are used only if fitted geometry 
is desired. There are 10 cards,) 

Format (5E 14. 8) 


upstream area coefficients. I = 1, 2, 3 
downstream area coefficients. 1 = 1, 2, 3 


(These cards are used only if standard nozzle 
geometry is desired. There are 2 cards, ) 
Standard geometries are wedge, cone or 
hyperbolic, axisymmetric nozzles and have 

area distributions of the form A = A + A ? ^+ 

2 1 L 
A % upstream of the throat, and 
^ 2 

A = Bj + B^oc + B^ downstream where % is 
nondimensional (x is measured from the throat 
and is negative upstream). Use of the fitted 
geometry permits the specification of seven 
different polynomials, two upstream of the 
throat and five downstream. The transition 
points must also be specified. For example 


x < ATP(l), A = Aj + A ^ 

ATP(l) £ x, ± 0, A = A 4 +A 5 x +A 6 *c 2 + 

3 4 5 6 

A 7 * + Ag* + A g* + A 1q *, 

0 4. * 4. ATP(2), A = Bj + B 2 * + B 3 * 2 

AT P(2) AT P(3), A = B 4 +B 5 * + Bg* 2 + 

3 4 5 6 

Byx + BgT + B 9 « + B 10 * 


ATP(3 ) 4 ATP(4), A=B n +B 12 * 

+ B , 7< 6 


! 
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ATP(4) 4 ATP(5), A=B j8 +B 19 ^ 


ATP(5) 


+ b 24 * 6 


A ~ B 25 + B zb* 


+ B 



Quadratic fits must be used for the first and 
third segments, since the program requires an 
explicit solution for ye in these regions. How- 
ever, these segments may be made very small, 
one or two integration steps. If a fitted geometry 
is used, all 10 upstream and 3 1 downstream 
coefficients must be read in. If fewer segments 
or lower order polynomials are required, 
zeros should be read in for the additional 
coefficients. The area distributions should be 
chosen to match as closely as possible at the 
transfer points. Slight discontinuities can be 
tolerated as long as the area distributions 
remains monotonic in both converging and 
diverging sections. (Also see SL on Card Type 4) 

Typical nozzle expansions span a wide range of temperature. This 
presents difficulty in accurately prescribing thermodynamic data. Simple- 
harmonic -oscillator formulas which are accurate at low temperatures may 
not be sufficiently precise at the higher temperatures. Polynomial fits of 
tabulated data (thermo-fits) may be used at higher temperatures. Consider- 
able flexibility has been built into the program and several options are 
available . 

1. All harmonic oscillator (HO) data. Set CTMXX = 0, NFIT = 0, 

KHO = 1. IGJ(J) need not be read in (Card Type 14). Only 
harmonic oscillator data need be read in. 
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2. All thermo-fit (TF) data. Set CTMXX = 0, NFIT = 1, KHO = 0, 
IGJ(J) = 1 for all species. Only thermo-fit data need be read in. 

3. Harmonic oscillator data for some species, thermo-fit and 
harmonic oscillator data for the rest, with a switch from TF to 
HO data at T = CTMXX. Set KHO = 1, ICJ(J) = 0 for HO species, 
IGJ(J) = 1 for the others. Set CTMXX equal to the desired 
switching temperature, which must be the same for all species. 

Set NFIT - 1, Only HO data need be read in for the HO species, 
while both sets of data are required for the rest. If the TF data 
is to be used for the entire temperature range, set CTMXX = 0 
and do not read in HO data for these species. 

It is essential that the two sets of data agree very closely at the switching 
temperature, otherwise oscillations may occur in the nonequilibrium solution. 


Card Type 14 Format (3011) 

IGJ(J) 1 if thermo-fit data is used for species J, 

J = 1 P . . . , ISS 

0 otherwise 

This card is used only if thermo-fit data are 
employed. 


The reading sequence for cards 15-18 is as follows: All thermodynamic 
data for species 1 is read in; Cards 15, 16, 17, 18 as required. Next, data 
for species 2, etc. Species are numbered in the order they appear in the 
oL j . matrix. 

Card T ype I 5 Format (4E15.7) 


TFA(J) 

a i 

TFB(J) 


TFC(J) 

c i 

TFD(J) 

d i 


(There will be 


IGJ =1.) 
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1 card for each species for which 



Format (3E15.7, A15) 


TFE(J) 


TFK(J) 


SHJAP(J) 


HP( J) 

species symbol 


(There will be 1 card for each species for 
which IG J = 1 . ) 

Card Type 17 Format (A6, 4E14.8,IZ) 


HP( J) 

species symbol 

ETAJ(J) 

V 

SBJ(J) 


THEVP(J) 

e' 


1 l 

SHJAP(J) 

A « 

IGM(J) 

number of electronic levels 


(There will be 1 card for each species. If 
CTMXX = 0 there will be one card for each 


species for which IGJ - 0. ) 

Card T ype 18 

Format (4(F3,0, E15#8)) 

CELJ (L, J), 
ELJ(L, J) 

If no electronic excitation is considered, one 
card must still be included, with GELJ (1, J) 


ground state degeneracy, ELJ(1, J) = 0, 

For each card of type 17, there will be enough 
cards of type 18 to describe the electronic 
levels, 4 to a card, a maximum of 10 levels 
permitted. 
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6. 2 Sample Calculation 


The sample case chosen to illustrate the use of the nozzle flow program 
is an expansion of ionized air through a contoured nozzle. It has been 
established that at the selected reservoir conditions, 7^ = 6500°K and 

- 1000 atm, the dominant source of electrons is NO . Thus NO~*~ is the 
only ion included in the calculation. An 8 species, 11 reaction chemical 
model is used with thermodynamic data in both harmonic -oscillator and 
thermo-fit form. The description of the necessary input data is summarized 
below and the corresponding input data cards are shown in Fig. 3. Sample 
pages of the output data are shown in the ensuing section to illustrate the 
printout formats. 

Summary of Input Specifications for Sample Calculation 

Run No. 1000 Upstream start, vibrational equilibrium 

Do frozen, equilibrium, and nonequilibrium 
solutions 

T,' - 6500 °K = 1000 atm. 

Species: e , N^, O^i A, N, O, NO, NO + 

Elements: N, O, A, e 

Reactions : 

k. 

1. N 2 +M^=2N+M M = N 2 

OxlO 21 ^' 1,5 exp |^-2. 2499 x 10 5 /R o T'J 

-A. 

2. N 2 +M^=:2N+M m=n 

A ( = 1. 5 X 10 22 T' " 1 * 5 exp [-2. 2499 X 10 5 /R o T'] 

3. N 2 + M 2N + M M = 0 2 , A, O, NO 

= 9.9 x 10 20 T' ' L5 exp [-2. 2499 x 10 5 /R o T / ] 
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A 

4. + M 20 + M M = °2 

Jk t =3.6x 10 2l T /_1,5 exp [-1. 1796 x 10 5 /R o T^| 

5. 0 2 + M 20 + M M = O 

= 2. 1 x 10 18 t' "• 5 exp [-1. 1796 x 10 5 /R o T^] 

4 . 

6. O z + M 20 + M M = N 2 , A, N, NO 

*4 = 1. 2 x 10 21 T /_1, 5 exp 1. 1796 x 10 5 /R o t] 

7. NO + M^=^N + O + M M = N 2> 0 2> A, N, O, NO 

Jk* =5.2x 10 Z1 T r1,5 exp [- 1 . 4996 x 10 5 /R o T^J 

8. 0 2 + N^O + NO 

Jk f =10 12 T / ' 5 exp [-6. 2 X 10 3 /R o r] 

A, 

9. N 2 + O ^NO + N 

Jki - 5 x 10 13 exp [-7. 552 x 10 4 /R o t] 

A* 

10. n 2 +o 2 ^no+no 

Jkf = 9. I X 10 24 T /_2 - 5 exp [-1. 2912 X 10 5 /r o t^ 

+ - k* 

11. NO + e ;= N + O 

Jk f = 1. 8 x 10 Z1 T , “ 1 * 5 

Nozzle Geometry = 1 cm, 

2 

Upstream: A = 1 + oL 

2 

Downstream: 0 ^ x <. 1. 25, A = 1 + . 6l41903x 

1.25 3.25, A = .6000286 + .63 9954 lx + 

.3582086* 2 
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3.25 C* ^ 5.08, A = 20.39065 - 15.40933*,+ 

4. 613686* 2 - .3664322* 3 
5. 08 <£ A = .006603241 + 1.607688* + 

. 3041264* 2 - .02642175* 3 + . 000948788* 4 - 

.001619982* 5 + .00001 0773 02sc 6 

Since the upstream fitted geometry involves only one polynomial the 
same coefficients are read in for both of the available upstream regions. 

An arbitrary transfer point, ^ = - 1 , is selected. Only 4 of the 5 downstream 
regions are utilized. The last transfer point is chosen beyond the range of 
the calculations. 

Thermodynamic Data 

Both thermo -fit and harmonic -oscillator data are to be used for each 
species, with the switch from the former to the latter description taking 
place at T' = 5000°K. The numerical values of the constants in the thermo- 
dynamic data may be seen from the input data cards in Fig. 3 or, in more 
compact form, in the program listing of the input data given in Fig. 4 in 
Section 6.3. The constants in the harmonic -oscillator data were taken from 
Reference 2. The coefficients in the thermofit data were taken from 
Reference 23. It has been established that the two thermodynamic 
descriptions are nearly the same at 5000°K. References 2 and 23 also 
contain the necessary harmonic -oscillator and thermo-fit data for the other 
ionized species present in high-temperature air (N + , 0 + , N 4 0 9 + , Ar + , O"). 

u Lf 

Constants and Test Values 

BZERO = 0 (no real gas correction) 

Maximum value of 4 , = 3 20 

Minimum value of temperature = 3 250°K (The solution was stopped at 
this point only for the purposes of the example. Usually the frozen-flow and 
equilibrium flow solutions are stopped at .01 of Tj ) 
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Initial An =0.01 
DELT2 = . 01 
DELT3 = . 001 


TTEST = . 05 
GTEST = . 10 
Interval Size Control = 4 


6. 3 Output Formats 

In this section the formats for printing the results of the machine 
computations are described. The results of the sample calculation, for 
which the input data were described in Section 0. 2, are used to illustrate 
the printout procedure. Then the diagnostic messages which may occur in 
the printout are listed and are used to indicate some of the experience which 
has been gained in using the program. 

For each computational case or run, the program lists the input data. 

The program listing of the input data for the sample calculation is shown in 
Figure 4. The format for this printout is evident from the figure. The 
specified values of the indicators IUPD and INEQV are denoted by printing 
UPSTREAM RUN or DOWNSTREAM RUN and VIBRATION EQUIL. or 
FROZEN VIBRATION on the input data listing. The program listing of the 
input data also includes the format for displaying the frozen, equilibrium 
and nonequilibrium solutions. If a given solution is not to be computed then 
the format for that solution is not printed. 

The results of the reservoir computation are printed after the program 
listing of the input data. The results of the reservoir computation for the 
sample air case are shown in Figure 5. The reservoir temperature, pressure, 
and density are in units of °K, atm, and gm/cc, respectively. The enthalpy 
is normalized by ^( 0 /R 0 T 0 > and the entropy is in cal/gm°K. The species 
concentrations are in moles /gm. 

The frozen flow solution, if computed, is printed after the reservoir 
values using the format given in Figure 4. The first page of the frozen-flow 
solution for the sample case is shown in Figure 6. The gasdynamic variables 
are all printed in dimensionless form. The entropy is always given in 
cal/gm° K. Notice that the printout of the composition is not repeated at each 
step in the frozen flow solution. 
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The throat conditions for frozen and equilibrium flow are shown in 
Figures 7 and 8 for the sample calculation* The first page of the equilibrium 
solution for the sample case is shown in Figure 9. The printout of the density- 
fit constants used in the starting procedure for the sample nonequilibrium 
solution are shown in Figure 10. If the solution is to be started in the 
downstream region then the conditions at the initial computational point 
downstream of the throat would be printed in place of the density-fit constants. 
The format for printing the starting values for a "downstream run" is the 
same as that for printing the equilibrium throat conditions. 

The printing procedure for the nonequilibrium flow results is 
illustrated in Figure 11. The first three pages of the printout for the 
sample calculation are shown in the figure and the format is that given in 
Figure 4. Again all of the flow properties are given in dimensionless form 
with the exception of the entropy which is given in units of cal/gm°K. The 
species concentrations are in moles/gm of mixture. In the nonequilibrium 
solution the current value of the indicator INEQ, which denotes when the 
numerical integration has begun, is printed at each step. 

The values of the nonequilibrium flow quantities which are printed for 
the steps where INEQ = 0 are actually the equilibrium flow variables plus the 
appropriate perturbations (see Figure 11). Also, when INEQ = 0, the 
equilibrium-flow temperature and species concentrations are printed at each 
step according to the format presented with the input-data listing in 
Figure 4. The above printout procedure permits the values of the perturba- 
tion quantities calculated in the starting procedure to be examined. 

After the first step on which the value of INEQ changes from 0 to 1 the 
values of the nonequilibrium flow properties printed are those obtained from 
numerical integration. The starting values for the numerical integration are 
the quantities on the step where INEQ is first equal to unity. Notice that 
after INEQ becomes 1 the printout of the nonequilibrium solution no longer 
follows the format shown in Figure 4, in that the equilibrium-flow temperature 
and species concentrations are not printed at each point. 
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All of the above output data are written on the system output-tape 
unit. For the nonequilibrium solution additional information is written on an 
auxiliary tape unit. The present program uses a tape unit designated logical 
tape unit 2, A page of these additional output data which was generated in 
the sample calculation is shown in Figure 12. For the steps corresponding 

to values of INEQ - 0 the quantities Xjl ' P z ' S h. ’ (bk ’ z = 1 * 2, 

. . . , /t , are printed. When INEQ = 0, the quantities ^ , P>. , and 

are evaluated using the equilibrium values of f ■ , T , p , andfr.) 

is evaluated using f. = /. + <5 and T»T + ST. At each step 

A v 7 o 

the run number and equilibrium -flow temperature are included in order to 

\ 

relate this additional printout to the printout of the flow properties and 
composition. The format for printing each step is as follows 


Run No, T 



s i 

s q 

i. «,)„ 

p > K 

X 5 6 

1. K CQ K 



i « 

• • 


The quantity X . evaluated using the infinite -rate equilibrium values 
. T , and p gives an indication of how closely the results of the 
Newton-Raphson procedure satisfy the condition of chemical equilibrium. 
As pointed out in Section 3.4 a comparison of (i an( 3 S t . is a 

check on the accuracy of the linearized perturbation calculation. 


The additional output which corresponds to steps where INEQ = 1 
contains the quantities r , , p. , p. x w 

V4, ** <C * yfc 

course, ^ > P- and are the values evaluated in performing the 

numerical integration of the nonequilibrium solution. Since there may be 
more reactions than species in a chemical model the fourth entry in each 
group of 4 quantities is meaningful only for the first s groups. The 
format for printing these quantities when INEQ = 1 is 


and 


dr i' d * 


Here, of 
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Run No . T 


i, p, Pl i, J* t t, £ p,i x t, p > it J*i 

t, ? 

*, p . p .t, ^ ^ 1 , — 




These quantities are useful in checking the progress of the non- 
equilibrium solution and in determining which reactions in the kinetic model 
are the dominant ones. 

In Section 5 . 3 an option for including the effects of moderate gas 
imperfections on the nozzle flow solutions was described. If this option is 
used and if the indicator IROBAR is read in as unity then the effective 
density, p , defined in Section 5.3 is printed with the output data. In the 
frozen and equilibrium -flow solutions p is printed in place of the entropy, 
which does not change in either case. In the nonequilibrium solution the 
correction is applied over the near -equilibrium portion of the flow and /& 
is printed in place of the Mach number. Notice that *p need not be printed 
even though this option is exercised. 

Diagnostic messages may appear interspersed with the output data 

obtained from the systems output tape. The diagnostic messages are of 

two types. The first includes error messages which are written by the 

24 

FORTRAN IV library system; the second includes diagnostic messages 
written by the program. 

24 

The error messages generated by the FORTRAN IV system are 
typically encountered when attempting operations such as taking the square 
root or logarithm of a negative number, taking the logarithm of zero, or 
taking the antilogarithm of a number greater than 88. The most frequent 
sources of these error messages have been eliminated as discussed below. 

In the nonequilibrium calculation a Mach number is computed from 


Equation (3 8). As can be seen from this equation M varies from 1 to oo as 
dJUhjd*/ din, p/d* varies from 0 to -1 . Since M is sensitive to small 

variations in this ratio of derivatives, an oscillation in the solution which 
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led to a value of ^dJM,p/d^< ^dl^A/dic^ would in turn result in taking 
the square root of a negative number. In the program this situation is 
avoided by calculating M from 


M = 


1.0 


/ j I +(cLJkv Afdt^ dl^p/ dt)\ 


The occurrence of this condition will be evidenced in the output by 
an oscillation in the Mach number. 


The program has been modified to eliminate FORTRAN IV error 

message printout in the calculation of the quantities £ # and \ i 

If one of the species concentrations on the reactant side of a reaction 

approaches zero and the reverse rate of a reaction becomes much larger 

^88 

than the forward rate, then operations such as jifuiO) and e are 

attempted. This situation is avoided for the current computational step by 
computing the reverse reaction-rate constant from the forward rate constant 
and the equilibrium constant, and changing the direction in which the 
reaction is written. 

The messages which may appear in the second type of diagnostic 
printout, i. e. that written by the program, are listed below. The 
significance and probable cause of each of the messages are described. 


(i) DISCRIMINANT NEGATIVE IN MAIN PROGRAM - This diagnostic 
refers to the calculation of from the quadratic formula (Eq. (44)), The 
sign used in the quadratic formula has been chosen based on the expected 
behavior of A (*) and hence this diagnostic message is seldom encountered. 

If encountered this diagnostic probably indicates an error in the input data. 

(ii) TOO MANY NEWTON-RAPHSON ITERATIONS - This diagnostic 
indicates that more than the allowed number of iterations have been performed 
in the Newton-Raphson procedure for calculation of the equilibrium com- 
position and pressure. It is obvious from the output data whether the failure 
occurred in computing the reservoir or a point in the equilibrium flow. 

Apart from inconsistencies in the input data this failure can occur when the 
various species concentrations have widely different values. One possible 
remedy would be to relax the test on AX^ which,as described in Section 2. 1, 
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defines convergence. See also the discussion given under subroutine 
SIMSOL in Section 7. 1. 


(iii) IN NE WRAP, CAPX(K) = 0, IN NEWRAP, P - 0. - These 
diagnostics are printed when the correction to the guess made for one of the 
unknowns in the Newton-Raphson iteration procedure would yield a value of 
zero for the next guess. This diagnostic usually occurs when the iteration 
procedure diverges and hence is usually accompanied by (ii). 


(iv) BETA MATRIX OF INSUFFICIENT RANK - As discussed in 
Section 4. Z, the matrix must be of rank (s- c) in order for the per- 

turbation quantities used in the starting procedure to be uniquely determined. 
Before performing the nonequilibrium flow calculation the rank of A 
checked and if it is not ( 5-c ) the computation is stopped. 


7 


is 


(v) DATEST AND DBTEST INCONSISTENT - This diagnostic refers 
to a failure in joining a nonequilibrium solution started in the upstream 
region to the desired solution in the downstream region. The failure referre 
to is that the test (DBTEST) on the percentage difference between the value 
of cCJh,A/ cLy calculated from the upstream region solution and that 
calculated from the specified A(ic) for the downstream region cannot be 
satisfied for the specified switching point (DATEST). If this diagnostic is 
encountered without (vi) then there is considerable nonequilibrium effect on 
the flow density upstream of the nozzle throat. This diagnostic usually 
indicates incorrect specified input data. See (vi)-(b) below - 


(vi) DLOGR IS POSITIVE - As discussed in Section 3. 3, a branch point 
exists slightly downstream of the geometric throat in the solution for a non- 
equilibrium nozzle flow. The primary purpose of inserting this diagnostic 
into the program was to indicate whether an upstream solution had been 
joined to the subsonic branch of a downstream solution. However, testing 
for points in the calculation where becomes positive has additional 

diagnostic value. The probable causes of failing this test at various points 
in the nozzle flow calculation are listed below. 



(a) Downstream run - If this diagnostic occurs in a calculation 
which is started downstream of the nozzle throat there are two 
possible causes. First, if the diagnostic occurs at a point away 
from the nozzle throat, then either an oscillation in the solution or 
decrease in the nozzle-cross - sectional area is responsible for 
dLL^p/d* becoming positive. An oscillation in the solution is 
accompanied by small integration step sizes. A discontinuity in 
dA/df, will usually occur at a transfer point in the A fa) 
relation. 

If this diagnostic occurs for a value of A very near 1 
(say i- 1. 005) then the numerical integration has been started too 
close to the nozzle throat. This situation is easily remedied by 
either halving the specified values of CECHII(I) and restarting the 
calculation as an upstream run or by doubling the values of 
CECHII(I) and restarting the downstream run, 

(b) Upstream run - Again there are two probable causes of this 
diagnostic. The first of these, which occurs away from the nozzle 
throat, is the same as the first type of difficulty listed above for 
the downstream run. If the diagnostic is encountered very near 
the throat and the numerical integration has begun, then the 
switch from the upstream region to the downstream region has 
been attempted too close to the nozzle throat. This situation is 
remedied by moving this switching point downstream, i.e. by 
increasing DATEST, (DATEST is initialized in Subroutine (INIT. ) 

(vii) TEMPERATURE GREATER THAN RESERVOIR VALUE - At each 
point in the nonequilibrium solution the value of the flow temperature is 
compared with the reservoir value. The test is intended as an over-all 
check on the progress of the calculation. This diagnostic indicates a 
radically incorrect computational result, most probably caused by an error 
in input data. 
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(viii) DISCRIMINANT NEGATIVE IN AXFIT (see i). 

(ix) MATRIX SINGULAR or INDEXING OR STORAGE FAILURE - 
These diagnostics are encountered in subroutine MATINV (see Section 7. 1), 
which is used to invert theoc matrix. The first diagnostic will occur if the 
rank of oc^ is n °t equal to c , i. e. if * s improperly specified. 

The second diagnostic will occur if the dimension statements of BTA and 
ALPIJ are not consistent, and will only be encountered when recompiling 
the program. 
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7. SUBROUTINES AND VARIABLES USED IN FORTRAN- 
LANGUAGE NOZZLE FLOW PROGRAM 

In this section the function performed by each subroutine in the pro- 
gram is explained. Flow charts are provided for the more complicated 
subroutines. Also, a list of all the FORTRAN variables is given, indicating 
their definitions and the subroutines in which they are used. 

7. 1 FORTRAN Subroutines and Their Functions 

Subroutines for which flow charts appear in Section 7. 2 are marked 
with an asterisk. After each description, the calling subroutines are 
written in parentheses. 

MAIN PROGRAM ''’ The routine which directs the computation is designated 
as the main program. This routine contains the logic for selecting 
which of the reservoir, frozen-flow, equilibrium-flow or nonequilibrium- 
flow computations are to be done on a given run. 

READ This subroutine performs the reading in of the input data. 

(MAIN PROGRAM). 

LIST This subroutine writes the input data on the system output tape. 

The format in which the input data are listed is given in Section 6. 

(READ). 

INIT Performs initialization and nondimensionalization of data. 

(MAIN PROGRAM). 

INTA First the specification of the chemical composition of the gas model 
is rewritten in terms of the designated independent species. Then the 
computation of the reservoir composition and thermodynamic properties 
is performed. (INIT). 

THERM This subroutine computes the thermodynamic properties, Jroj , 

for each species at a given temperature. 

RAP, NONEQ, PERT, RNKT). 


/v ’ , and 

(INTA, FROZEN, NEW 


0 

V 
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FROZEN This is the subroutine which controls the calculation of the 
throat conditions for frozen flow and the frozen-flow solution, 

(MAIN PROGRAM). 

PROP At each computational step in the frozen-flow calculation this sub- 
routine calculates the pressure from Equation (25), the enthalpy from 
Equations (12) and (84b), the density from Equation (lid) and the 
velocity from Equation (15). (FROZEN). 

PRTA This subroutine controls the printing of the output data for the 
nonequilibrium flow solution. The output format is discussed in 
Section 6. (NONEQ). 

NRMAX This subroutine controls the computation of the equilibrium-flow 
throat conditions, (MAIN PROGRAM). 

NEWRAP The Newton-Raphson iteration procedure for finding the 

equilibrium-flow conditions at a given temperature, entropy, and 
reservoir condition is contained in this subroutine. (NRMAX, 

EQUIL, NONEQ, AXFIT). 

EQUIL This subroutine controls the calculation of the equilibrium-flow 
solution (MAIN PROGRAM). 

NONEQ This is the controlling subroutine for the nonequilibrium 

solution. It also contains the logic for the starting procedure and for 
switching from an upstream to downstream region. (MAIN PROGRAM). 

COMM This subroutine computes the production terms in the species 

conservation equations (Eq. (2 9)). If the nonequilibrium integration 
has been started the pressure and entropy are also computed using 
Eqs. (lib) and (16), respectively. (NONEQ, RNKT). 

GEOM In this subroutine the area ratio at a given nozzle station is 

computed from the specified polynomial relation. If the nonequilibrium 
integration is started upstream of the nozzle throat the density is 
computed from the density-fit relation and the area ratio is obtained 
from the continuity equation. (COMM). 
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* df' aLT 

EXACT In this subroutine the derivatives of --l b- and are computed 

for use either in the numerical integration of the nonequilibrium 
equations or for use in the perturbation calculations. (NONEQ, RNKT). 

PERT This subroutine computes the perturbation quantities used in 

the starting procedure for the numerical integration of the nonequilibrium 
equations . (NONEQ). 

AXFIT The function of this subroutine is to invert the A(d)r elation to 

obtain the value of ^ corresponding to a given area ratio in the 
starting procedure of the nonequilibrium flow calculation. (NONEQ). 

THROAT When the nonequilibrium solution is started upstream of the 

nozzle throat it is obtained as a function of the equilibrium flow 
density. This subroutine modifies the specified geometry to match 
that computed for the nonequilibrium flow upstream of the nozzle 
throat. (NONEQ). 

RNKT ' This subroutine performs the numerical integration of the nonequili- 
brium flow solution. (NONEQ), 

MATINV This subroutine transposes the o£ x y matrix for the purpose of 
specifying the chemical composition of the mixture in terms of the 
independent species. The calling sequence for this subroutine is 
CALL MATINV (BTA, ISC, 64) where BTA is the matrix to be trans- 
posed, ISC is the size of the square matrix, and 64 is the 2nd dimension 
of this matrix according to the dimension statement. (INTA). 

SIMSOL In both the equilibrium and nonequilibrium flow calculations it Is 
necessary to solve a set of linear, simultaneous algebraic equations. 

This subroutine is one which is available for such a function in the 
program library at CAL. This subroutine is written in MAP language 
and hence is peculiar to the IBM 7044 version machine in use at CAL, 
i.e. one including double precision hardware. The calling sequence 
for the subroutine used in the nozzle flow program is CALL SIMSOL 
(AA, ISSP2, 22). The array denoted by AA consists of the matrix 
of coefficients and the right hand sides of the simultaneous equations 
to be solved. The right hand sides of the equations are listed in the 
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last column of AA. The second number, here ISSP2, in the calling 
sequence denotes the size of the matrix of coefficients. The third 
number, here 22, is the number of rows allowed for in the dimension 
statement for AA. 

One of the computations in the program which uses this sub- 
routine is the Newton-Raphson procedure for finding the composition 
at a point in an equilibrium flow. In an equilbrium nozzle expansion, 
the concentrations of some species, for example atomic species, 
become quite small. As discussed in the section on input data, species 
which are expected to vanish in an equilibrium expansion to low tempera- 
ture should not be chosen as independent species. In early work with 
the program the Newton-Raphson iteration was found to diverge even 
when the concentrations of the dependent species became vanishingly 
small. This situation was remedied by using a SIMSOL subroutine 
which employs double precision arithmetic. While this accuracy 
is not consistent with the single precision arithmetic used throughout 
the program this procedure permits continuing the computation even 
when dependent species concentrations underflow the digital limit of 
the machine and hence are set to zero. This subroutine included 
in the program is one using double precision arithmetic. If such a 
subroutine cannot be used, some provision must be made to set to 
zero the concentrations of dependent species when they become less 
than a given small number. Otherwise the equilibrium solution cannot 
be continued to temperatures where species become vanishingly small. 

As discussed in Section 5.1, in considering ionized flows in 

equilibrium, the electrons, which are an independent species, may 

vanish. The electron concentration cannot be allowed to underflow 

and hence in this case the gas model is truncated and all ionized species 

are dropped from the calculation when the electron concentration 

-30 

becomes less than 10 . When the concentration of the electrons 

- 15 

becomes less than about 10 moles per gram the accuracy of the 
results for the ionized species is questionable. If desired the 
electrons could be dropped from the calculation at a larger concentration. 

(INTA, NEWRAP, EXACT, PERT). 
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SAVE Frequently it is necessary to terminate program execution for 
considerations of time and/or to check the progress of the solution. 

This can be done by storing the contents of the computer memory on 
magnetic tape. Since the subroutine used for this purpose is peculiar 
to the CAL system, a dummy subroutine has been inserted in the 
deck. Users may wish to insert their own subroutine for this 
purpose. (NONEQ). 

SSWTCH, PDUMP, EXIT These are standard subroutines available in 

. “ ' 2 5 

systems utilizing ^FORTRAN IV language programming. 

7.2 Flow Charts of FORTRAN Subroutines 

This section contains flow charts of the more complicated subroutines 
in the nozzle flow program. These flow charts do not depict every statement 
in the FORTRAN Source deck but rather connect the description in the text to 
the actual program. The following symbols are used in these charts: 

Operational block. This symbol refers to 
computational instructions. When appropriate 
the description in an operational block contains 
references to equations in the text. The num- 
bers above the upper right-hand corner refer to 
the numbers of the cards in the FORTRAN source 
deck to which the block refers. 

Decision block. This symbol denotes points 
in the program beyond which two or more paths 
are possible. The information in the block 
questions the value of some indicator. Again 
the number outside the block refers to the card 
in the FORTRAN program. 

Calling block. This symbol indicates a call 
statement for a FORTRAN subroutine. 
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Connector. This symbol is used to display a 
connection for a case where arrows cannot be 
used conveniently. The number in the circle 
refers to the number of the card in the FORTRAN 
deck corresponding to the instruction to which 
the connection is made. 
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MAIN PROGRAM 


430 
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SUBROUTINE INTA 


3820-3840 








SUBROUTINE THERM 


5030-5050 


FORM VARIABLES NEEDED 
IN COMPUTATION OF SPECIES 
THERMODYNAMIC PROPERTIES 


IS IGJ(J) = 0? 


IS CT < CTMXX? 


5080 

\yes 


5180-5530 


51 10-5160 


COMPUTE -fij , Af , AND CAUSING 
THERMO-FIT DATA IN EQUATIONS (79) 
AND (80) 




IS INEQV = 0 ? 







SUBROUTINE FROZEN 



77 










SUBROUTINE NRMAX 



78 











SUBROUTINE NEWRAP 



79 












SUBROUTINE EQUIL 



80 





















SUBROUTINE CO M 



82 









IS X >ATP( I) ? 



SET HI = 2. N2 = 3. AND SI = B, 


1290 0-12920 


SET HI = 5. N2 ' 10. AND SI = 8 U L_. 


12950.12970 


SET Nl = 12. N2 = 17, AND SI * Bi 


13000- 13020 


SET Nl - 19. N2 = 2N, AND SI = B, 


12790-12870 


N2 


COMPUTE 


SI FROM SI = SI + Xjl B jX 

N2 


+ 1-N1 


AND S2 FROM S2 =Z ( i-l) 8 X 1 
Nl ' 


i -Nl 


12880 


IS IUPO = 0 ? 


YES 


1 3060- I 3 MO 


13120 


dA d£ nA 

, AND FROM SI ANO S2 

d X d X 

IMPUTE p FROM EQUATION (13) 


RETURN 
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SUBROUTINE EXACT 



84 







SUBROUTINE PERT 



85 
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SUBROUTINE THROAT 



87 










SUBROUTINE RNKT 


15890-15950 













YES 


1 6620-16630 



CALL THERM 
CALL COMM 
CAL! FX ACT 



168 5 0 - 168 7 0 



RETURN 



















7 . 3 FORTRAN Variables Used in Nozzle Flow Program 

The variables which appear in the COMMON statement of the nozzle 
flow program are listed in this section. The FORTRAN symbols, their 
definitions, and the subroutines in which they are used are given. The 
numbers used to denote the subroutines in which a FORTRAN variable is 
used refer to the following list. Where possible the corresponding algebraic 
symbol is used to define a FORTRAN variable. 



1 . 

MAIN 

12. 

EQUIL 


2, 

READ 

13. 

NONEQ 


3. 

LIST 

14. 

COMM 


4. 

INIT 

1 5. 

GEOM 


5. 

INTA 

16. 

EXACT 


6. 

THERM 

17. 

PERT 


7. 

FROZEN 

18. 

AXFIT 


8. 

PROP 

19. 

THROAT 


9. 

PR TA 

20. 

RNKT 


10. 

NRMAX 

21. 

MATINV 


1 1 . 

NEWRAP 



FORTRAN Var 

iable Subroutine 

Definition 

A 


1, 15 


o c , constant in density-fit 
relation 

AA(I, J) 


5, 1 1, 17 


matrix element 

AAA(I, J) 


16 


matrix element 

ACOM(I) 


2 


comment which appears in 
list of input data 

AFNTS 


1, 13, 18 


used to denote area ratio 


calculated for equilibrium 
flow in starting procedure 
for nonequilibrium solution 


AF NX 7, 9, 1Z, 13, 15, 19 A » area ratio ( A= a!/A* ) 
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FORTRAN Variable 

Subroutine 

Definition 

AL.PIJ(I, J) 

2, 3, 5, 16, 17 

oc^j, number of atoms of j in 
element per molecule of the 
i*h species 

AMACH 

7, 9, 12, 16 

hi , Mach number 

AR 

13, 19 

indicator which denotes where 
becomes positive and 
number of times switch from 
upstream to downstream region 
has been attempted 

ARBA 

4, 13 

maximum number of tries at 
switching from upstream to 
downstream region 

ARBB 

4, 13 

counter for number of times 
ups tr earn - downs tr earn 
switching point is moved 
downs tream 

ASUB(I) 

1-3, 15, 18 

A, , coefficients in area-ratio 
expression for upstream region 

ATP(I) 

2, 3, 15, 18, 19 

transfer points in area-ratio 
expre s sion 

B(I, J) 

13 

dummy variable used in 
ranking jg , . 

BCHI(I) 

9, 17 

1, evaluated for £ + 57'>T+&T>/5 + 
he / / 

BE(I) 

13, 17 

A ■ A A, 

BET(I) 

5, 11, 16 

t) : -1 where >*• =5~ jA . 

TV 

BETA(I, J) 

13, 14, 16, 17 

Pif ’ A- ; - = ~ 

BLBK(I) 

13 

dummy variables for down- 
stream area coefficients 

BSUB(I) 

1-3, 13, 15, 18, 19 

, coefficients in area- 
ratio expression in downstream 
region 

BTA(I, J) 

5 

dummy variable used in 

transposing u- • 

*a 
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FORTRAN Variable 

Subroutine 

Definition 

BZERO 

2, 3, 5, 8, 1 1 

^ see Eq- ( llc ) 

C 

1,15 

C , constant used in density- 
fit relation 

CAI(I) 

2, 3, 13, 14 

A- j constant factor in 
(Eq. (82)) 

CAPQ(K) 

2, 3, 5 

Q , number of gram-atoms of 
* element k 

CAPX(J) 

1, 5, 10, 11, 18 

X - » mole fraction of j** 1 species 

3 

CAPXTH(J) 

1, 10, 18 

vT , mole fraction at throat 

CARB 

13 

dummy variable for A(j£) 

CCI(I) 

5 

intermediate variable for 
computing molecular weights 
of species from molecular 
weights of elements 

CCPJ(J) 

6, 14, 17 

’ AA'j Mr 

CDIJ(I, J) 

5, 1 1, 16 

t) ** , stoichiometric coefficients 
in equilibrium formation 
reactions 

CEACT(I) 

2,3,14 

E-act- 9 activation energy 

CECHII(I) 

2,3,13 

6^ ♦ » test on which 

determines where numerical 
integration is begun 

CGI(I) 

5, 11 

molecular weights of species 

CGMU(I) 

5, 11 

JU Xi 

CH 

1, 5,8, 11, 14 

/-/ , enthalpy 

CHA 

1, 5, 8, 1 1, 14 

, reservoir or stagnation 
enthalpy 

CHI(I) 

9, 14 

^ , defined by Eq. (31) 

CHII(I) 

5, 11 

intermediate variable in 
computation of equilibrium 
constant based on mole -fractions 

CLNIMC(I) 

14 

M'-o 
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FORTRAN Variable 


Subroutine 


Definition 


CLNPI(I) 

14 

JUP. 

CLNT 

6, 14 

2r»T' 

CM 

5, 11, 14 

7^ , molecular weight 

CMA 

1, 5, 7, 8, 1 1, 13, 14, 16, 17 

7fl e .reservoir molecular 
weight 

CMW(I) 

2,3,5 

7f[~ , atomic weights of 

elements 

CRA 

4, 5, 14 

R 0 , universal gas constant 
(=1. 98647 cal/mole °K) 

CRP 

4, 13, 14, 17 


CRRB 

5, 7, 11 

intermediate variable in 
calculation of entropy 

CRS 

13 

dummy variable for entropy 

CSTA 

13, 14 

°- 5 *~(K/ R X) 

CT 

1, 4, 6-14, 16-18, 20 

T , tempe rature (t = J */Tj^ 

CTAP 

1-6 

_ / 

T 0 f reservoir temperature 
in °K 

CTB 

13, 20 

dummy variable for T 

CTC 

17 

dummy variable for T 

CTMAX 

1, 7, 10, 18 

_ X 

T , throat temperature 

CTMXX 

2-4, 6 

1 

temperature for switching 
from thermo-fit to harmonic- 
oscillator model for species 
properties 

CTP 

6 

T / , temperature in °K 

CTPL 

4 

JU,Tj 

CTT 

9, 13 

temperature at previous step 
in nonequilibrium calculation 

CX 

1,9, 13, 15, 18-20 

, distance along streamtube 
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FORTRAN Variable 

Subroutine 

Definition 

CXB 

13, 20 

dummy variable for & 

CXMAX 

2, 3, 13 

maximum value of ^ desired 
in nonequilibrium solution 

DATEST 

4, 13 

value of A at which non- 
equilibrium solution is 
switched from upstream to 
downstream region 

DBTEST 

4, 19 

percentage difference allowed 
between calculated and specified 
dLJLnjfK/d^t a t switching point from 
upstream to downstream region 

DELT1 

1, 4, 7, 10, 12, 13 

temperature increment used 
in frozen flow and equilibrium 
flow calculations 

DELT2 

1-3, 18 

temperature increment used 
to start nonequilibrium solution 
in the downstream region 

DEL.T3 

2, 3 

temperature increment used 
to start nonequilibrium solution 
in upstream region 

DELTAX 

2, 3, 13, 19, 20 

Ad 

DGJ(J) 

9, 16, 17, 20 

dlj /d*. 

DLOGA 

13, 15, 16, 19 

dAvA/oid 

DLOGR 

13,15, 16 

ctjlrv p /d& 

DT 

16, 20 

dT/d* 

ELJ(JL, J) 

2-4, 6 

t'Jtj > energy of£^ electronic 
level of j th species 

ELMENT(K) 

2, 3 

chemical symbols for elements 
in mixture 

ENT 

7,8' 

intermediate variable used in 
frozen flow calculation 

ETAI(I) 

2, 3, 14 

, temperature exponent 

in A - 
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FORTRAN Variable 

Subroutine 

Definition 

ETAJ(J) 

2-4, 6 

/J . number of atoms in j 
specie s 

FLUX 

7,8,10-13,18 

/° u 

GELJ(L, J) 

2, 3, 6 

fy. degeneracy of elec troni 

level of jth species 

GJ (J) 

1,9-14, 16, 17, 20 

species concentration in 
moles / gm 

GJA(J) 

1, 5, 7, 8, 10, 12 

(lj) 0 species concentration 

in reservoir 

GJB(J) 

13, 20 

dummy variable for 

GJC(J) 

17 

dummy variable for 

GTEST 

2, 3, 13 

allowed value of At*/?* in 
Runge-Kutta integration step 

HDELX 

13,19, 20 


HP(J) 

1-3 

chemical symbol for species 
in mixture 

IC 

2, 3, 5, 11 

number of ions included in 
chemical model 

IGJ(J) 

1-3, 6 

indicates whether thermo-fit 
data is to be used for the jth 
specie s 

IGM(J) 

2-4, 6 

number of electronic levels 
for jth species 

IM 

1, 4, 11 

index which is incremented 
when electrons are dropped 
from calculation 

INEQ 

9, 13-16 

indicates when numerical 
integration of nonequilibrium 
solution has begun 

INEQV 

2, 3, 6 

indicator for selecting whether 
equilibrium or frozen vibra- 
tional model is to be used 
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FORTRAN Variable 


Subroutine 


Definition 


IP 

7,9, 12, 13 

printout line-counter 

IROBAR 

2, 3, 7, 9, 12 

indicates that effective density 
is to be printed out 

IRUN 

1-3, 9 

run number 

ISC 

2-5, 10, 1 1, 13, 16, 17 

C , number of elements in 
mixture 

ISC PI 

13, 16,17 

C + 1 

ISMC 

1, 4, 5, 13 

5-C 

ISMCNR 

1,4, 11 

current value of fs-c) in 
equilibrium flow calculation 
(decreased by IC when elec- 
trons are dropped from 
calculation) 

ISR 

2, 3, 9, 13, 14, 16, 17 

A* , number of reactions 

ISS 

1-14, 16-18, 20 

5 , number of species in 

mixture 

ISSNR 

1,4, 11 

current value of s in equilib- 
rium flow calculation 
(decreased by IC when elec- 
trons are dropped from 
calculation) 

ISS PI 

13, 16, 17, 20 

5 ♦ 1 

ISSP2 

13, 16, 17 

s + Z 

ISSP3 

13, 16, 17 

s * 3 

ISSP4 

13,16 

s + 4 

ISW1A 

1 -3 

indicates whether frozen flow 
solution is desired 

ISW2A 

1 -3 

indicates whether nonequilibrium 
flow solution is desired 

ISW3A 

1-3 

indicates whether equilibrium 


flow solution is desired 
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FORTRAN Variable 


Subroutine 


Definition 


ISW4A 

1.2 

indicates whether another 
calculation is to be done 
after completion of one case 

ISW6A 

1.2 

indicates when only reservoir 
calculation is to be done 

ISW5A, ISW1B 
ISW2B, ISW3B 
ISW4B, ISW5B 
ISW6B 

2 

additional indicators not 
currently used in program 

ITB{I) 

13 

dummy variable for indicators 

IUPD 

1-3, 13, 15, 16 

indicates whether calculation 
is in upstream or downstream 
region 

IZERO 

1 -3, 5, 6, 13 

logical variable used to denote 
zero 

JJK 

1.4, 11 

indicates when electrons are 
dropped in equilibrium calculation 

KHO 

2-4 

indicates whether any harmonic- 
oscillator data are used 

KKUR 

2 

indicates whether the third 
body matrix is used in the 
chemical model 

KUR(I, J) 

2, 3, 14 

. , third body matrix 

“tf 

LC 

13, 18 

indicates the first time sub- 
routine AXFIT is called 

Ml 

5, 10, 1 1 

index equal to c *+ t 

NAFIT 

2, 3, 15, 18, 19 

indicator for selecting 
standard or fitted geometry 

NFIT 

2-4 

indicates whether any thermo- 
fit data is to be used 

NIT 

1 1 

counter for number of iterations 
in Newton-Raphson procedure 

NNN 

13, 20 

integration step counter 
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FORTRAN Variable 

Subroutine 

Definition 

NNS 

13,20 

number of times A has been 

increased by current value of 
SC 

NQS 

2, 3, 13 

number of successful steps 
before A is increased 

NTEST 

4, 5, 1 1 

number of iterations allowed 
in Newton-Raphson procedure 

PCT 

9, 13, 17 

ST , perturbation in tempera- 
ture 

PC TEST 

4, 13 

tolerance on £ when testing 
size of 6^ . ^ 

PERTGJ(J) 

9, 13, 17 

6/- > perturbation in species 
concentrations 

PGJ(J) 

1 1 

intermediate variable in 
Newton-Raphson calculation 
of equilibrium composition 

PI(I) 

9, 13, 14, 17 

P, defined by Eq. (32) 

PICHI(I) 

9, 14, 16 


PRES 

1, 7-12, 14, 18 

•fy, pressure * ■p'/p'*') 

PRESA 

1-3, 5, 14 

/ 

, reservoir pressure in 

atm 

PRESB 

7, 12 

dummy variable for 

PRESTH 

1,7, 10, 18 

it 

pressure of throat 

PRHO 

17 

5 f> perturbation in density 

QM(I) 

5, 11 

number of moles of j th 
components in one mole of 
mixture 

QQ(I) 

2, 3, 14 


RHAP 

1, 5, 14, 17 

pj reservoir density (gms/cc) 

RHO 

1 , 7-15, 17 

P , density (/° - p'/pf} 

RHOB 

7, 12, 13 

dummy variable for p 
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FORTRAN Variable 


Definition 


Subroutine 


RHOBAR 

7-9, 11, 12 

RHOC 

17 

RHOP 

14 

RHPL 

14 

RHTH 

1, 7, 10, 15 

ROBARA 

5, 7, 8, 11, 12, 14 

ROBARP 

5,8, 11 

SAJ(J) 

4, 6 

SBJ(J) 

2-4 

SC 

13, 20 

SCPG 

14, 16 

SDCHI(I) 

9, 13, 17 

SDELTX 

13, 19 

SDGJ(J) 

13, 20 

SDT 

13, 20 

SEN 

1,5, 7, 9, 12-14 

SENT(J) 

5, 6, 8, 1 1, 14 

SHDELX 

13, 19 

SHJ(J) 

5, 6, 8, 1 1, 14, 16, 17 

SHJA(J) 

4,6 

SHJAP(J) 

2-4 


, effective density defined 
by Eq. (84) 

dummy variable for p 
^'density in gms/cc 

s&7V p' 

p density at throat 

H'M' 

^'effective density in reservoir 

(tj constant in harmonic - 
oscillator expression for yU*. 

Jrj constant defined by Eq.(76) 

factor by which integration 
interval is changed 

s u perturbation in 

initial value of A 

A fj , change of in 
integration interval 

AT , change of T in inte- 
gration interval 

, entropy in cal/gm °K 

, species entropy, 5^*- 3y 

initial value of A v/2 

Ai , species enthalpy 

J / * 

X , formation enthalpy 

, formation enthalpy in 
cal/ mole 
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FORTRAN Variable 


Subroutine 


Definition 


SHPG 

8, 14 

r i 

SKIL(I) 

5, 11 

intermediate variable in 
calculation of dependent 
species concentrations from 
concentrations of components 

SL 

2,3,13,18 

ji characteristic length 

SL64 

13, 14 

l /6,470 

SM 

1,7, 10, 12, 13, 15, 18 

H critical mass flow rate 

SS 

7 


SU 

1 , 7-9, 11,12, 14, 15 

U , velocity u = a '/%/R o T 0 ' 

SU2 

8,11, 14, 16, 17 

U* 

SUMG 

1 4, 16 


TB( J) 

13 

dummy variable for quantities 
necessary to return to previous 
step in nonequilibrium calcu- 
lation 

TEMPI 

1,18 

discriminant in calculation of 
from quadratic formula 

TEST 

4, 5, 11 

size of correction to guesses 
used to define convergence of 
Newton-Raphson procedure 

TESTB 

4, 7, 10 

tolerance allowed in locating 
temperature at throat for 
frozen and equilibrium flow 

TFA(J) 

2, 3, 6 

CLj , coefficient in thermo- 
fit expression for JL^ and /^J/T 

TFB(J) 

2-4, 6 

Jlrj y coefficient in thermo-fit 
expression for Jhj and ^uj/T 

TFC(J) 

2-4, 6 

Cj > coefficient in thermo-fit 
expression for Ji . and/y/r 

TFD(J) 

2-4, 6 

dj > coefficient in thermo-fit 
expression for and yttj/T 

TFE(J) 

2-4, 6 

ej, f coefficient in thermo-fit 
expression for J^j and 


99 



FORTRAN Variable 

Subroutine 

Definition 

TFK(J) 

2, 3, 6 

J ^ coefficient in thermo-fit 
expression for Jt,- and /T 

THEV(J) 

4, 6 

© v . characteristic vibrational 

temperature /T c ') 

THEVP(J) 

2-4 

characteristic vibrational 
te&ipe rature in °K 

TPRINT 

2, 3, 9, 13 

temperature interval at which 
results of nonequilibrium 
calculation are to be printed 

TSTOP 

2-4, 7,12,13 

minimum value of temperature 
desired in solution 

TTEST 

2, 3, 13 

allowed value of AT/T in 
Runge-Kutta integration step 

UP 

13,18 

indicates when nonequilibrium 
solution which has been started 
upstream is restarted down- 
stream of nozzle throat 

XMJAT(J) 

5, 6, 1 1, 14, 17 

Of 0 

Mj/T where jul • is species 
chemical potential^! 

XNUI(I) 

13,14 

i 1 

stoichiometric coefficient 
of jth species on reactant side 
of itb reaction 

XNUIJ(I, J) 

2, 3, 13, 14 

XNUIJP(I, J) 

2, 3, 13 

stoichiometric coefficient 
of species on products side 

of i^h reaction 

ZP 

8, 1 1 


ZPA 

5, 7, 11 
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APPENDICES 


MODIFICATIONS OF THE BASIC NOZZLE FLOW PROGRAM 
FOR ALTERNATIVE APPLICATIONS 


The basic computer program calculates the expansion from an equilibrium 
reservoir through a specified converging-diverging nozzle. Modifications of 
the program to extend its applicability are described in this section. The 
changes in the governing equations and in the specification of the input data 
are given. As for the basic program, the copies of the FORTRAN source 
cards are available for these modified versions. The modifications were 
carried out with an emphasis on convenience rather than programming 
efficiency. Consequently these modified versions contain unnecessary 
variables and logic which could be removed if desired. 
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APPENDIX A 


NONEQUILIBRIUM EXPANSION FROM ZERO VELOCITY 
THROUGH A SPECIFIED PRESSURE DISTRIBUTION 

In the first modification of the program the pressure distribution rather 
than the streamtube geometry is specified. The program is then capable of 
computing the flow along a streamline for a given pressure variation. Typical 
applications of this program include the design of a nozzle with a specified 
pressure distribution and the computation of the nonequilibrium flow around 
the surface of a blunt body with an assumed surface pressure distribution. 

Since the frozen and equilibrium flows are independent of the rate of 
expansion they are computed in the same way as when the nozzle geometry is 
specified. The properties of a frozen or equilibrium flow at a given point 
are then determined by the pressure at that point. For programming conven- 
ience the procedure for locating the maximum in f>u. when the streamtube 
cross section is specified is not altered. The area ratio calculated at each 
point in the frozen or equilibrium flow then gives the spreading of the stream- 
tube relative to its minimum cross-section. 

In Section 3. 1, the governing equations for a nonequilibrium expansion 
through a converging-diverging nozzle were presented. In the computer 
program these equations are reduced to (5+2) first-order differential 
equations, or in finite -difference form to (stZ ) algebraic equations for the 
derivatives cLlj/di.* dT/dv > and djhvpfdx* The equations for the slopes are 
in terms of the flow properties at a given station and the value of dJb, A / * dlt 
for the specified geometry. For the expansion through a specified pressure 
distribution these equations are written in terms of d'fj/di, $ dT^d# > di*» hjd4 > 
and dJinp/dit . 

The equations for the conservation of elements ( Z 8 ) and conservation of 
species (30) are written in the same form for a specified pressure distribution. 
In addition, dtl/df, is eliminated from the momentum equation using the energy 
equation. Also introducing the equation for the enthalpy^ (12), yields 


104 



( 86 ) 






dT 

d$. 



d *p. ; 


ctx. 


Then, using the energy equation together with (12) to eliminate the velocity 
in the continuity equation and combining the result with the logarithmically- 
differentiated state equation yields 


( y r )dr ,^-(ud7n s u % dUA u* 

[\T jn / fTi \ % **)**'% d* ' % d* 

Equation (87) may be simplified using Equation (86). The result is 


y , _L *L 

fr, ** V ** 


/ dtvvA. 


l f, 122. \ dlrup; 

9q\ 7r(u.y d* 


(87) 


( 88 ) 


which together with Equations (28), (30), and (86) comprises the necessary 
set of [s+2 ) equations. 


Once the flow properties and composition are known at a given point, 

the derivatives and can be calculated. Then these values are 

<t* at* 

numerically integrated to obtain the temperature and composition at a new 
station. Having determined T and the s the density may be obtained from 
the state equation using these values and the specified pressure at the new^ , 
Also, the enthalpy may be computed from Equation (12), the velocity from 
Equation (15), and, using the equilibrium critical mass flow, the streamtube 
cross-sectional area is found from Equation (13). The area ratio of the 
nonequilibrium streamtube flow will not be unity at the point of minimum 
cross-section. However, this calculation again is only intended to indicate 
the relative spreading of the streamtube and is not a necessary part of the 
computation. 


When the nonequilibrium nozzle-flow equations are written in terms of 
a specified pressure distribution there is no branch point associated with the 
solution. Consequently, it is not necessary to modify the present computation 
upstream of the point where the velocity is equal to • Thus, 

could be eliminated from Equations (86) and (88) since there is no good reason 
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for carrying it along (as there was for retaining in the equations for a 

specified geometry). However, retaining in Equation ( 88 ) greatly 

simplifies the modification of the program to take a specified pressure 
distribution. Also, computing ■ at each point permits a simple calculation 

for a Mach number based on the "speed of sound, a~ — • 


M = 



dlxsK/ctlC 

d Jh'p'/ d-# 


(89) 


In addition to the Mach number the entropy is also computed at each step using 
Equation (16). 

The nonequilibrium calculation for a specified pressure variation is 
started in the same manner as that for a nozzle flow. The nonequilibrium 
solution is considered to be a perturbation about the infinite - rate equilibrium 
solution and the numerical integration is started when the perturbation 
calculation indicates significant departure from equilibrium. 


The perturbation calculation must be modified for the case of a specified 
pressure distribution. The perturbed form of the element conservation 
equations remains valid. 

Y~ oC A <5 i- = 0 'k. - I, Z, • • • C (55) 

; r ' 3 

In order to obtain the perturbed form of the species conservation equations 
the £ ’ 5 should be written as ^ s ‘ Then 


d7- 


where again it is assumed 
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(90) 


Since is prescribed, S-p>= 0. Evaluating the other partial derivatives 

in (90) leads to 
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Since the density is not necessarily very near the equilibrium value, 
the condition <5S = 0 rather than S F = 0 is used to replace the perturbed 
momentum equation. From the definition of the Gibbs free energy 

f'= h'-t's ' 

the perturbation in F to first order is seen to be 


8 F = &H - $ ST -T (s S') 


(92) 


6 9 


L *i 


Also, since F r F(j,p>, , then also to first order^ 

SF = W) 5T + |r) s r+tit) 

r > V r X r f /= ' '} i. 

or, using 0 this equation becomes 

8 F = -SST + 7?l.t. /*■*&** 

f-i • • 

Now combining (92), (93), and the first-order perturbation of (1Z) gives, in 
the case of 6 S' = 0, 


(93) 




(94) 


dJUA 


In rewriting the nonequilibrium flow equations the term ^ ^ — was 
retained to simplify reprogramming. The perturbation in the streamtube 
cross-section is retained here for the same reason. Using the perturbed 
forms of the continuity, energy, and state equations together with 5i» = 0 
leads to 

<SA ^ (KZ; . -'h., . (% • 


T-ZFP-Wt?} SS-t) 

Equation (94) may be used to simplify Equation (95) to 

<5 A X. M,\ „ . ST 


8 T 


(95) 




(96) 
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In the program the equation 

5 , = hist*- = s r _ 


T 


5. /y 


(97) 


is employed to write jul: in terms of 5 * and - . Then Equations ( 9 a) and 

/ $ o 

( 9 h) become 


I T . <5T =0 

X. - T (sS -Mv (y 7ri ' 


44- = + + 

A T £7 \ ( u L7 



(98) 


(99) 


Now Equations (55), (91), (98), and (99) are a set of ( s + 2 ) algebraic equations 
fo r 67 ^. , 6 T , and £ A . 

The procedure described in Section 3.Z for starting the nozzle-flow 
computation is directly applicable here. The equilibrium solution is obtained 
at successive temperature steps from the reservoir. At each point the 
equilibrium slopes are computed and used to determine the perturbations in 
T and the . When the perturbations are sufficiently large the numerical 

integration of the nonequilibrium solution is begun. In the computation of the 
equilibrium slopes the ( s-c ) equations obtained by differentiating the con- 
dition of chemical equilibrium (see Section 3.2) should be written in terms 


of 


d7 x. 


d.T 


c£ A 


and 


d Jtn, 


* d'i 1 — dbtf — 5 ctllva — * This can be accomplished by using 

the differentiated form of the continuity equation and the momentum equation 
whereupon Equation (85) becomes 

y~ 1 _ Qi—ll — 

fe 1 k d* " A* T* Al T7 d* T d * 

_/o 1) , Q±zh ' _^L£l -n 

V* i ' d# pu d'jc ( 100 ) 

In the starting procedure for the nozzle flow computation the value of 
A(€) for the equilibrium solution is used to find the corresponding /£ from 
the specified area distribution. In the starting procedure for the streamtube 
flow the equilibrium flow pressure is used to find jC from the specified 
pressure distribution. The method employed to invert PM is the same 
as that described in Section 3. 2 for inverting AU). 
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As mentioned above, the modification of the computer program to 
permit a specified pressure distribution was carried out with a minimum of 
revision. As a result, the specification of the input data for this case is 
identical to that described in Section 6 with the following exceptions. The 
absence of a singularity at the throat of the streamtube excludes the "upstream” 
region distinction made in the nozzle flow case. Consequently the logic in the 
specified-pre s sure -distribution program is that for the "downstream” nozzle- 
flow region. Accordingly, the indicator IUPD is always read in as 0 for this 
case. Furthermore, the coefficients in the specified -fvfa) relation are read in 
using the locations for the area- relation coefficients, but starting with the 
first downstream region. That is, *6 is measured from the stagnation point 
rather than from the nozzle throat. The storage locations available for 
specifying the upstream area distribution have been used to provide an additional 
interval for specifying a pressure distribution. Thus, the polynomial relations 
given in Section 6 for a fitted geometry become for a specified pressure 
distribution 


Ip'U) - 8 , + B z 4 + B 3 4 % 

fW -- £ B. / 

Jl- U-n ) 

i»(*) = YL 

i -a 

2 4 y . . 

-?(*) * Z_ B x * 

1x18 


0 £ « 4 ATP(2) 
ATP(2) 4ATP(3) 

ATP(3) 4 . ATP(4) 

ATP(4) <=■ H 4. ATP(5) 

ATP(5) £ 4 4 ATP(l) 

* i ATP(l) 


In the starting procedure the quadratic formula is used to invert in the 

first interval, 0 4^ 4 ATP(Z). Since behave s differently than A(t) the 

opposite sign must be used in the quadratic formula when inverting % ). 

Thus, in this case Equation (44a) must be written as 
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(44b) 


in the program. 


/fi 


2 . 6 3 


The program changes to effect the calculation for a streamtube flow 
with a specified pressure distribution involved the following changes in the 
FORTRAN variables listed in Section 7. 


AFNTS > calculated from equilibrium solution 

DLOGA ^ GijUsf,/ ci # 

DLOGR > dd^A/d* 

The subroutines in this modified version of the program are the same 
as described in Section 7 except for those in which the changes discussed in 
this section are made. The changes in these subroutines are briefly 
summarized below. 

MAIN PROGRAM The computation of the first point in the nonequilibrium 
solution uses the logic for starting the calculation at the reservoir. The 
option of starting the nozzle flow computation at the throat is eliminated. 

Also, the density fit computation is omitted. 

NONEQ All of the logic associated with switching the calculation from an 
upstream region to a downstream region is eliminated. 

COMM The computation of the pressure from the state equation is eliminated. 

GEOM The computation of A and dl^Ajcli, for the intervals in the upstream 

region are eliminated. The computation of the density from the density-fit 
relation is also eliminated. The computation of A and dd^k/dt is modified to 
calculate ^fx and dJLrbjxfd'l . The computation of the density from the state 
equation is added. 
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EXACT The changes discussed above for the governing equations for the 
derivatives d^/dd , dT/dy, and ^/«/^^are made. The changes incorporate 
the FORTRAN notation change given above. The formula for the Mach number 
computation is changed to that given in Equation (89). 

PERT The modifications of the perturbation computation to account for a 
specified pressure distribution are effected in accordance with Equations (91), 
(98), and (99). The perturbed state equation is used to relate 8p to ST and 
the Sfj's. The computation for the in terms of , ST , and *P is 

then kept the same as in the nozzle flow calculation. 

AXFIT The logic for determining the proper set of coefficients for inverting 
-pO#) is changed to account for the difference in the behavior of and A(d) . 

Since there are no upstream and downstream regions in the nonequilibrium 
calculation for a specified pressure distribution, subroutine THROAT is 
eliminated. 
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APPENDIX B 


NONEQUILIBRIUM FLOW FROM A NONEQUILIBRIUM 
INITIAL CONDITION THROUGH EITHER A SPECIFIED 
PRESSURE OR AREA DISTRIBUTION 


The basic nozzle-flow program and the modification described in 
Appendix A both deal with flows which start from an equilibrium stagnation 
region. The modification described here is designed to compute a quasi- 
one -dimensional, nonequilibrium flow from a finite -velocity, nonequilibrium, 
initial condition. Either a given pressure distribution or a given streamtube 
area distribution may be the specified boundary condition. This version of 
the program is applicable, for example, to streamtubes in a body flow field 
for which the pressure distribution is known from either measurements or an 
equilibrium flow-field calculation. Also, the chemical- relaxation zone behind 
a strong shock wave can be solved by specifying the conditions behind a 
translational- rotational- vibrational equilibrium shock and a constant stream- 
tube cross-section. Furthermore, this modification can be used to continue 
a calculation originally started from a stagnation point with either boundary 
condition. 

When the nonequilibrium solution is started from a nonequilibrium initial 
condition the starting procedure discussed in Section 3.2 becomes unnecessary. 
If the flow properties and composition are given at an initial point the derivatives 
may be computed and the numerical integration begun immediately. The 
governing equations for the nonequilibrium flow through a streamtube of 
given cross-section are Equations (28), (30), (34), and (36) of Section 3.2, 

As discussed in Appendix A, if the streamtube pressure distribution is 
specified Equations (86) and (88) replace (34) and (36). The numerical 
integration of these equations proceeds as discussed in Section 3. 1 for the 
specified area distribution and in Appendix A for the specified pressure 
distribution. 

The modified version of the computer program described here was also 
obtained with a minimum of reprogramming. Consequently, the input format 
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is again essentially the same as described in Section 6 for the nozzle-flow 
program. The values of T 0 and p)/ to be specified are respectively the 
temperature and pressure at the given initial condition. In addition to the 
input data described in Section 6 the local values of all the species concentra- 
tions in units of moles per gram, the local temperature nondime ns ionalized 
by T c , the nondimensional axial distance, the total enthalpy in cal/gm. and 
the molecular weight at the given initial condition must be specified. It is 
important when using this version of the program to continue a calculation 
that the original values of T c , and '-p'J , and the co r re sponding molecular 

weight are specified. Otherwise the computed values in the continued solu- 
tion will be normalized by different factors from those used in the first part 
of the solution. 

This version of the program uses the logic for the downstream region 
of the basic program. Thus, the indicator IUPD should be read in as zero. 
Five intervals are available for specifying A(x) and ; if more are 

needed they can be obtained by repeated starting of this version of the pro- 
gram. One of the available indicators, ISW1B, is used in the program to 
select either the governing equations for a specified A(x) or a specified 
If A(<t) is specified this indicator should be read in as zero. If •p'fcjis 
specified it should be read in as unity. When -f>0z) is specified the indicator 
ISW1B is also used to effect the changes in the FORTRAN variable notation 
listed in Appendix A. 

The number of gram-atoms and the molecular weight of each element 
need not be given because all the concentrations are specified at the initial 
point. These cards should be eliminated from the input data deck. Also, 
the input cards for the upstream area coefficients should not be included. 

The variables and tests which must be specified for the frozen and equilib- 
rium calculations may be left blank on the appropriate cards (see Section t>). 

When this version of the program is applied to computing the flow 
behind a shock wave, for either a specified A(x) or , the initial value 

of A Ti' should be reduced to 10 -6 cm. Also, the integration tests, GTEST 
and TTEST should be relaxed to about 0. 2 and 0. 1, respectively. 
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Due to the elimination of the computations for the reservoir, frozen 
flow, and equilibrium flow, and the procedure to starting in equilibrium the 
following subroutines are eliminated in this version of the program; 

INTA NRMAX PERT 

FROZEN NEWRAP AXFIT 

PROP EQUIP MATINV 

Also, the logic from the nonequilibrium calculation for a downstream region 
is removed from NONEQ and incorporated in a modified MAIN PROGRAM 
peculiar to the present version. Since the logic associated with starting 
from equilibrium and with switching from an upstream to a downstream 
region are unnecessary, the remainder of subroutine NONEQ and subroutine 
THROAT are eliminated. This modification of the program also involved 
changes in the following subroutines; 

READ The instructions for reading in the number of gram-atoms and 
molecular weights of the elements and the A(%) or 00 coefficients for an 
upstream region are eliminated. The instructions for reading in the initial 
species concentrations , temperature, axial location, and H 0 are added. 

The procedure for preparing the additional input data is as follows; 

First the Y - ^ are read in according to the format used on card 
$ 

type 13 of Section 6 . Since there may be as many as ZO species there may 

be 4 cards to read in the initial . The values of hfj (in cal/ gm) , T } 

$ 

and In are read in after the 7 . ^ , again using card type 13. As mentioned 
0 * 

above, when using this version of the program to continue a solution, the 
specified value of should be the molecular weight of the mixture at the 

conditions and T o 

LIST Changes appropriate to the differences in input data are made. 

COMM The changes indicated for subroutine COMM in Appendix A are 
included as an additional option when is prescribed. 

GEQM The computation A(V) or -p'fc) is allowed for and the logic for 
an upstream region is eliminated. 
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EXACT The computation of the derivatives needed in the nonequilibrium 
solution for either a specified Afa) and are included. If (V.) is 

prescribed, Equations (86) and (88) are used rather than (34) and (36). The 
calculation of the equilibrium values of the derivatives is eliminated. The 
Mach number is computed from the formula appropriate to either a specified 
A fa) (Eq. (38)) or a specified pte) (Eq. (89)). 
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APPENDIX C 


NONEQUILIBRIUM EXPANSION FROM AN EQUILIBRIUM, 
FINITE -VELOCITY INITIAL CONDITION 


The modification of the program discussed in this section also treats 
the case where the initial flow velocity is nonzero. The problem treated 
here is the expansion of a flow assumed to be initially in equilibrium, but 
moving at supersonic speed, through a divergent channel. This version of 
the program is applicable to the flow in non-reflected type shock tunnels 
where the flow behind the incident shock is assumed quasi-steady. 

Exercising the option of including ionized species also allows the program 
to be applied to the flow downstream of an MHD accelerator. 

As for the nozzle expansion from an equilibrium reservoir state, the 
initial equilibrium state is found for a specified temperature and pressure. 
Since the initial state is not a stagnation point the total enthalpy of the flow 
must also be specified in order to evaluate the velocity. The frozen and 
equilibrium expansions from this initial point can also be obtained following 
the same procedure as discussed in Sections Z.Z and Z.3. In this case, 
however, the mass flow is known at the initial point. The cross-sectional 
area is referred to that at the initial point and the origin for is taken 
where A fa) = 1, i.e. the initial point. 

Since the computation using this program assumes the flow to be in 
equilibrium at a finite velocity, the nonequilibrium solution must be started 
differently than the nozzle expansion from a stagnation point. The 
perturbation starting procedure discussed in Sections 3.Z and 3.4 is valid only 
when the flow starts from an equilibrium stagnation point where the non- 
equilibrium flow actually approaches the infinite-rate equilibrium flow. 

The assumption made in this section that the flow is in equilibrium at a 
finite velocity is a physical idealization and the starting procedure is also 
artificial. The starting procedure adopted here is similar to that described 
in Reference 4 to start the numerical integration of a nonequilibrium nozzle 
flow. 
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Starting values for the integration for the nonequilibrium solution are 
obtained by assuming that over the first interval, A^ s * the species 
gradients are all zero. In other words, the flow is assumed frozen. Then 
at ^ r A* s the flow will be at some nonequilibrium state which approximates 
the actual flow. The numerical integration from then on proceeds exactly as 
discussed in Section 3. 1. 

The selection of an initial value of A% is influenced by the numerical 

difficulties of the type discussed in Section 3.5. In the present application 

the best initial value must be determined by trial and error. Experience 

with the use of this program for airflows indicates that an initial A of 
-3 

10 cm is a satisfactory value. 

The input data for the version of the program discussed herein is 
identical to that described in Section 6 for the nozzle flow calculation. In 
addition the total enthalpy of the flow must be specified in cal/gm°K. Since 
there is no nozzle throat and hence no upstream and downstream regions the 
indicator IUPD should be specified as 0. Although the coefficients of A>C#) 
in the upstream region are not used in this calculation, the appropriate number 
of blank cards must be included in the input data deck. Also, notice that 
ATP(2) is the first dividing point between the first and second intervals in 
the specification of A . 

As for the other two modifications of the basic nozzle-flow program 
the present version was obtained with emphasis on minimizing the amount 
of reprogramming. Consequently the subroutines in this program are the 
same as described in Section 7 for the basic program with the following 
necessary exceptions. 

MAIN PROGRAM The instructions for calculating the equilibrium throat 
and the density-fit constants are eliminated. 

READ The instruction for reading in the total enthalpy in cal/gm°K is 

added. h/ is read in after all the input data discussed in Section 6 and 
using the same format as card type 13, 
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INTA The computations of the velocity at the initial point and of the total 
mass flow are added. 

FROZEN The logic for locating the maocimum in and hence the throat 

for frozen flow is eliminated. 

EQUIL The procedure for starting the equilibrium calculation from a finite 
velocity initial state is incorporated. 

NONEQ The logic associated with the perturbation starting procedure and 
the switch from an upstream to a downstream region are eliminated. 

COMM The calculation of the rates for the individual reactions is 
eliminated on the first step, since <L fj / d,#, - 0 for this interval. The 

indicator LC is used to determine the first step. In the basic nozzle-flow 
program LC is used to indicate the first step in subroutine AXFIT, which is 
unnecessary in this program. 

GEOM The calculations of the f (#,) , A(&) , and a/cC /t for the 

upstream region are eliminated. 

The following subroutines are eliminated in this version of the 
program. 

NRMAX 

PERT 

AXFIT 

THROAT 
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N ATOM CONCENTRATION _ N ATOM CONCENTRATION (MOLES PER GRAM OF MIXTURE) 




Figure 2 COMPARISON OF RESULTS OBTAINED FOR NITROGEN - ATOM CONCENTRATION 

USING FOURTH-ORDER RUNGE-KUTTA AND MODIFIED RUNGE-KUTTA INTEGRATION 
PROCEDURES. 
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Figure 3a INPUT DATA CARDS FOR SAMPLE NOZZLE FLOW CALCULATION 
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Figure 3c INPUT DATA CARDS FOR SAMPLE NOZZLE FLOW CALCULATION 
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Figure 3f INPUT DATA CARDS FOR SAMPLE NOZZLE FLOW CALCULATION 
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Figure 3g INPUT DATA CARDS FOR SAMPLE NOZZLE FLOW CALCULATION 
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Figure 3h INPUT DATA CARDS FOR SAMPLE NOZZLE FLOW CALCULATION 
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Figure 3i INPUT DATA CARDS FOR SAMPLE NOZZLE FLOW CALCULATK 
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Figure Ma PRINTOUT OF INPUT DATA FOR SAMPLE CALCULATION 
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l.OOOOOOE 00 0. 1.00000 OE 00 l.OOOOOOE 00 0. l.OOOOOOE 00 0. 

0 . 0 . 


DOWNSTREAM AREA POLYNOMIAL COEFFICIENTS 


l.OOOOOOE 00 

0. 

0. 

0. 


0. 

0. 

6.6032*lE-03 

0. 


8.1*19036-01 
2.0390656 01 
1.6076886 00 


6.000286E-01 
-1.5*09136 01 
3.0*126*6-01 


6.3995*16-01 
*.6136866 00 
-2.6*21756-02 


1.5820866-01 

-3.66*3226-01 


0. 

0. 


9.4878806-0* -1.6199826-01 


0. 

0. 

1.077102E 05 


U> 


AREA FIT TRANSFER POINTS 


-l.OOOOOOE 00 1.250000E 00 1.250000E 00 5.080000E 00 


l.OOOOOOE 01 


SPECIES 

ATOMS 

PER MOLECULE 


CHEMICAL CONSTANT 

CHAR. VIBRATIONAL TEMP. 

ENTHALPY OF FORMATION 

ELEC TRONIC 

LEVELS 

fc- 


l.OOOOOOE 

00 


-1.A928J2E 01 


-0. 



-0. 

1 


M2 


2 . 0000006 

00 


-*.2163006-01 


1.3512*06 03 


-0. 

* 


02 


2 - OOOOOOE 

00 


1.07*5006-01 


2.218970E 01 


-0. 

5 


AR 


l.OOOOOOE 

00 


1 .8655706 00 


-0. 



— 0. 

3 


N 


l.OOOOOOE 

00 


2.986800E-01 


-0. 



1.1259066 05 

5 


0 


l.OOOOOOE 

00 


*•9120006-01 


-0. 



5.8980006 0* 

6 


NO 


2. OOOOOOE 

00 


5. 19*1006-01 


2.6991806 01 


2.1*77006 0* 

1 


NO* 


2. OOOOOOE 

00 


1.7861006-01 


1. 17295 OE 01 


2.1531006 05 

5 


SPECIES 

1 DEGENERACY, ELECTRONIC 

ENERGY LEVEL 1 








t- 

2. 

-0. 











N2 

1 . 

-0. 


1 . 

l .*168506 05 6. 

1.70*7506 

05 l. 

1 .715600E 

05 




02 

5. 

-0. 


2. 

2.261700E 0* 1. 

1. 7725006 

0* 3. 

1.0119806 

05 

3. i.*Z3900E 05 



AR 

1 . 

-0. 


5. 

2.6610706 05 1. 

2.680420E 

05 






•4 

* # 

-0. 


6. 

S.*97*00E 0* *. 

5.512500E 

0* 6. 

8.2*55006 

0* 

12. 2.3821006 05 



0 

5. 

-0. 


1 . 

*. 5*62006 02 l. 

6. *898006 

02 5. 

*.5168006 

0* 

1. 9.6615006 0* 

5. 2.109070E 

05 

NO 

*. 

-0. 


2. 

1 .257000E 05 *. 

1.1128106 

05 






NO* 

1. 

-0. 


6. 

1.1521206 05 1. 

1*5898706 

05 6. 

2.0873206 

05 

2. 2.095900E 05 




SPECIES 

A 


B 

C 

0 

E 

K 

H6AT Of FORMATION 

£ - 

2.5000006 

00 

-0. 

-0. 

-0. 

-0. 

-1.1735006 

oi -o. 

N2 

3. *51*836 

00 

3.0883126-0* 

-*.251*286-08 

2. 7392956-12 

-5**683206-17 

3.0712696 

00 -o. 


Figure 4b PRINTOUT OF INPUT DATA FOR SAMPLE CALCULATION 



132 


02 

*4 

0 

*0 

NO* 


3.24947JE 00 
2.S63282E 00 
3.006922E 00 
2.S94143E 00 
3.7562L66 00 
3.3973SSE 00 


4.963449E-04 
-3.591770E-05 
■J.UW25E-04 
-5. 00 891 4E -05 
2 • 09 396 1 E -04 
3.74939*6-04 


-6. 7017SJE-0B 
7 • 4692 09E - 09 
6. 3118136-08 
l . I 99502E-09 
-/.6JWIE-08 
-6.062030E-08 


4.4433396-1? 
-6.747034E-11 
- 4 . 1652036 - 1 ? 
- 8 . 6016116-13 
1 ■ 6903 32 E- 12 
4. 637506E- 1 2 


-1 . 00028 1 E- 1 6 
2.2340196-17 
9. 3 34896 E- I 7 
2.1 481008*1 7 
-J.611323E-17 
-I.I07704E-16 


5.915022E 00 
4« 000939E 00 
1.103476E 00 
4.600615E 00 
3. 61 1 16 7E 00 
4. 2005636 00 


- 0 . 

- 0 . 

1.1259066 05 
5. 8980006 04 
2.I4770OC 04 
2 . 353300C 05 


OUTPUT FORMAT FOR FROZEN FLOtf 

SOLUTION 




TEMPERATURE PRESSURE 

DENSITY 

VELOCITY 

MACH NUMBER 

ENTROPY 


AREA 


OUTPUT FORMAT FOR EQUILIBRIUM 

SOLUTION 




TEMPERATURE PRESSURE 

SPECIES CONCENTRATIONS 

OEMS 1 TV 

VELOCITY 

MACH NUMBER 

ENTROPY 


AREA 


OUTPUT FORMAT FOR NONfcQulL I BRIOM SOLUTION 

TEMPERATURE PRESSURE DENSITY VELOCITY 

TERPEAAf URE--SPEC l£S CONCEMTRATI ONS 
EQUILIBRIUM TEMPERATURE— EQUILIBRIUM VALUES 


MACH NUMBER 


ENTROPY 


AREA 


A INEQ 


Figure 4c PRINTOUT OF INPUT DATA FOR SAMPLE CALCULATION 



133 


RESEiwain-ioao 


TEMPERATURE* 6530*00 OE^SITr* 3.Q6905906 PRESSURE- 1000.000 

E flTHALPV" 5.111* EMTftUPV* 2.1602 MOLECULAR *6 1 - 2 6. 1 652 

£•" 2.292E-06 

HZ 2.*M£-0 2 
02 1.61M-01 

Aft 1.219E-0* 

H *.6ALfc-0* 

O 6. ftftSE-01 
MO 4.T*«6*01 
MO* 2.2926*06 


Figure 5 SAMPLE PRINTOUT FOR RESERVOIR CONDITIONS 



FROZEN FLOW Silt Ilf ION 


UJ 

4^ 


I.OOOOOE 00 
9.90000E-0L 

9. 800006- 01 
9.70000E-01 
9.60000E-01 
9.50000E-01 
9.60000E-01 
9.30000E-01 

9.2 ooooe-oi 

9.1 00006- 01 
9.00030E-01 

8.900006- 01 

8.800006- 01 
8.70000E-01 
8.600006-01 
l.)0090l>t>l 
8.600006-01 
8. >00006-01 

а. zooooe-oi 

8. 100006- 01 
8.00000E-0! 
7*900 OOC-Oi 
7*800006*01 
7*700006-01 

7. 800006- 01 

7.500006- 01 
7*800 006-01 
7.30000E-0L 
7*200006-01 

7.100006- 01 

7.000006- 01 

6.900006- 01 

б. 800006-01 

6.700006- 01 
6.600006-01 

6.500006- 01 
6.600006-01 
6.300006-01 
6.200006-01 

6.100006- 01 

6.000006- 01 

5. 900006- 01 

5.800006- 01 

5.700006- 01 
5.600006-01 

5.500006- 01 
5.600006-01 
5.30000E-01 
5.20000E-Q1 
5 « lOOOOE-Ol 
5.000006-01 


i.oooooe oo 

9.58201E-01 
9. I 7781E-01 
8. 78T05E-01 
8.609606-01 
8.066566-01 
7.692126-01 
7. 351866-01 
7.023396-01 
6.706656-01 
6.600736-01 
6.105956-01 
5.821806-01 
5.568026-01 
5.286316-01 
5.030626-01 
6. 786086 -0 l 
6.551026-01 
6.326996-01 
6.107736-01 
3.899016-01 
3.698586—01 
3.50620E-01 
3.321666-01 
2 . 966166-01 
2.80715E-01 
2.656586-01 
2.51226E-01 
2.37601E-01 
2.261656-01 
2.1 1501E -01 
1.99390E-0I 
i. 878176-01 
1.767666-01 
L. 662156-01 
1.561566-01 
1.665666-01 
1.376316-01 
1 .287396-01 
1*206766-01 
1.126196-01 
1 .051626-01 
9*808 7 3E -02 
9* 1 38206-02 
8.503236-02 
7.902506-02 
7.136706-02 
6 • 7985 76 -02 
6. 292876-02 
5.816336-02 
5 • 36 79 36 - 02 


I.OOOOOE 00 
9.67880E-01 
9*365116-01 
9.058826-01 
8.759806-01 
8*6679)6-01 
8.183106-01 
7.905216-01 
7.63612E-01 
7.369736-01 
7.111936-01 
6.860626-01 
6.615696-01 
6.377036-01 
6.166556-01 
5.918156-01 
5.697716-OL 
5. 683156-01 
5.276376-01 
5.07128E-01 
6. 873766-01 
6.68176E-01 
6.695126-01 
6.313816-01 
3*900216-01 
3*762876-01 
3.589976-01 
3.661656-01 
3. 297266-01 
3.157266-01 
1.021666-01 
2.889726-01 
2.762026-01 
2.638276-01 
2.518616-01 
2.602376-01 
2.29007E-01 
2. 181656-01 
2.076666-01 
1.976986-01 
1.876996-01 
1.782606-01 
1.691166-01 
1.603196-01 
1.518636-01 
1.636826-01 
1.358286-01 
1.282 75t-0l 
1. 210176-01 
1 . 1 6067 1 -0 k 
1.073596-01 


0. 

2.916936-01 
6.121566-01 
5.066926-01 
5.826636-01 
6.513216-01 
7.133616-01 
7.703856-01 
8.236366-01 
8.732606-01 
9.203236-01 
9.650866-01 
1.007866 00 
1.06882E 00 
1.088256 00 
1.126276 00 
1.163026 00 
1.19866E 00 
1.233206 00 
1.266806 00 
1.29952E 00 
1.331616 00 
1.362566 00 
1.392 966 00 
1.662 1 36 00 
1.670356 00 
1.698026 00 
1.525186 00 
1.55186E 00 
L. 578066 00 
L.60383E 00 
L.62918E 00 
1 . 656 L 36 00 
1.678706 00 
1.702906 00 
1.726756 00 
1.750276 00 
1.773656 00 
1.796336 00 
1.818916 00 
1.661206 30 
1.863206 00 
1.886966 00 
1.906626 00 
1.927ft56 00 
1.9686)6 00 
1.969376 00 
1.989886 00 
2.010176 00 
2.030266 00 
2.050106 00 


0. 

2. 561 728 8E -01 
3.66015S5E-01 
6. 67985866-01 
5.19826396-01 
5.86066866-01 
6.63025666-01 
6.98077 33E-01 
7.50119666-01 
7*99768996-01 
8*67680386-01 
8.93591856-01 
9.38378116-01 
9.82065036-01 
L. 02677166 00 
1.06669516 00 
1.10796906 00 
1.16863006 OO 
1.1 888301 E 00 
1.22862786 00 
1.2681036E 00 
l. 307 1002E 00 
1* 36629916 00 
1.3851)196 00 
1.556)9126 00 
1.6768)516 00 
1.5166579E 00 
1.5525)536 00 
1.590503)6 00 
1 .62859726 00 
1.66685026 00 
1.70529216 00 
1.76)95726 00 
1.78287826 00 
1.822082)6 00 
1.86160826 00 
1.90168106 00 
1*9617)866 00 
l *98261 036 00 
2.02352866 00 
2.06513316 00 
2.10725676 00 
2.1699306E 00 
2.19319966 00 
2.237L015E 00 
2.2816721E 00 
2.32696326 00 
2.37*01196 00 
2.61986566 00 
7.66757996 00 
2.51619516 00 


2.160176 

00 

2.16017E 

00 

2.160176 

00 

2.16Q17E 

00 

2.160176 

00 

2.16017E 

oo 

2.16017E 

00 

2.16017E 

00 

2.160176 

00 

2.16017E 

00 

2.160176 

00 

2.160176 

00 

2. 1601 7E 

00 

2.160176 

00 

2. 1601 76 

00 

2.160176 

00 

2.160176 

00 

2.160176 

00 

2*160176 

00 

2.160176 

00 

2.16017E 

00 

2.160176 

00 

2-16017E 

00 

2.160176 

00 

2.16017E 

00 

2.160176 

00 

2.160176 

00 

2.160176 

00 

2.160176 

00 

2. 1601 76 

00 

2.160176 

00 

2.160176 

00 

2.16U17E 

00 

2.160176 

00 

2* 160176 

00 

2.160176 

00 

2.160176 

00 

2.160176 

00 

2.160176 

00 

2.160176 

00 

2.1 60 1 7t 

00 

2.160176 

00 

2.160176 

00 

2.160176 

00 

2.160176 

00 

2.160176 

00 

2.160176 

00 

2.160176 

00 

Z. 160176 

00 

2.16017E 

00 

2.160176 

00 


0. 


2.371366 

00 

1.73)296 

oo 

1*663366 

00 

1.310796 

00 

1-213036 

00 

1.16609E 

00 

1.098566 

00 

l. 066286 

00 

1.039596 

00 

1.022166 

00 

1.010656 

00 

1.00)616 

00 

1.000286 

00 

1.000526 

00 

1 .003 7)E 

00 

1.0096 l E 

00 

t -017956 

00 

1.028596 

00 

1.06160E 

00 

1.05633E 

00 

1.073)16 

00 

1.0923)6 

00 

1.11339E 

00 

1.189666 

00 

1.215686 

00 

1.266066 

00 

1.276626 

00 

1.307516 

00 

1.362806 

00 

1.380616 

00 

1.621086 

00 

1 « 66636E 

00 

1.510616 

00 

1.560016 

00 

1*612786 

00 

1.669166 

00 

1.729)66 

00 

1.79)656 

00 

1.862606 

00 

1.9)5906 

00 

2.016556 

00 

2.098766 

00 

2.188976 

00 

2.285716 

00 

2.389566 

00 

2.501096 

00 

2.621056 

00 

2.750216 

00 

2.88966E 

00 

3.039716 

00 


Figure 6 SAMPLE PRINTOUT FOR FROZEN FLOW SOLUTION 
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■'W 


FROZEN THROAT* 1000 

TEMPERATURE- 0*045840 MASS FLOA« 0.4690 DENSITY- 3.627157 PRESSURE « 0.S4JM1 


Figure 7 SAMPLE PRINTOUT OF THROAT CONDITIONS FOR FROZEN FLOW 



TE«PEftArtffc€* 0.903994 HASS FLOH* 0.6S59 OENSITV- 0. 62*117 PRESSUHf* 0.S57679 
EHfHALPV • A. SSI! VELOCITY* l.OSOJE 00 


e- i. owe- 06 
N2 2.A5AE-02 
02 l * 820E-0) 
AR 3.219E-0A 
H 2.46U-04 
0 6. 2S0E-0) 

HO 4*S9*E-0J 
HO* 1.099E-06 


Figure 8 SAMPLE PRINTOUT OF THROAT CONDITIONS FOR EQUILIBRIUM FLOW 
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EQUILIBRIUM SOLUTION 


L.OQOOQE 00 L .OOOOOE 03 1 • 00000 E 00 3.000036-99 0.00000006- 19 2. 1 60 176 00 0.000006- 19 

*. 2926-06 2.4J3E-02 1.4156-03 3.2196-04 4. 8816-04 6.8886-01 4.7686-01 2.2926-06 

9.900306-01 9.619206-0 1 9.545236-01 5.388056-01 5.09956296-01 2.160176 00 2.028216 00 

2. 1586-06 2. 6566-02 1.6536-05 5*2196-06 6.5776-06 6.8596-03 6.7566-01 2.1986-06 

9.803306-01 8.906196-01 9.10640E-01 6.791176-01 6. 15 l 96 786 - 0 l 2.1601T£ 00 1.509156 00 

1.9926-06 2.6186-02 1.6866-09 9.2196-06 6.2856-06 6.7756-03 6.7986-09 1.9926-06 

9.700306-31 8. 991776-01 8.682996-01 5.867916-01 5.16016686-01 2.160176 00 1.287)66 00 

1.8696-06 2.6606-02 1.5266-09 9.2196-06 6.006E-06 6.715E-03 6.7216-01 1.8596-06 

9.600 006-01 7. 907116-01 8.27651E-01 6.775976-01 6.22602156-01 2.160176 00 1.169866 00 

1.721E-06 2.6626-02 1.5666 01 1.2196-06 1.7996-06 6.6596-01 6.7066-03 1.7216-06 

9.500306-01 7.449206-01 7. 08048E-01 7.576606-01 7.00911996-01 2.160176 DO 1.098586 00 

1.5956-06 2.6446-02 1. 6056-03 9.2196-06 9.6846-04 6.587E-01 6.6896-03 1.5956-06 

9.400306-01 7.001196-01 7.500446-01 S.90062E-01 7. 71915146-0 l 2.160176 00 1.059546 00 

1.6766-06 2.6466-02 1.6486-01 1.2196-06 9.2616-06 6.5196-03 6.4716-03 1.476E-06 

9.100306-01 6.580026-01 7.199976-01 8.96728E-01 8.19059526-01 2.160176 30 1.025116 30 

1.9696-06 2.4696-02 l. 6936-01 1.2196-06 9. 0096-06 6.6686-03 6.6526-01 1.369E-06 

9.200006-01 6.179046-01 6.780666-01 9.588596-01 9.02825606-01 2.160176 CO 1.008866 00 

1.2576-06 2.6516-02 1.7606-01 9.2196-06 2.7906-04 6.9746-01 6.6926-01 1.2576-06 

9. 100006-01 5.797996-01 6.440076-01 1.017906 00 9.69969606-01 2.160176 00 l. OOLITE 00 

1.1566-06 2.4596-02 1.709E-03 9.2196-06 2.5016-04 &. 2976-01 4.6UE-09 1.156E-06 

9.000006-01 5.494906-01 6.L1187E-01 1.0T267E 00 1.02902726 00 2.160176 00 1.00067E 00 

1.0626-06 2. 6556-02 1.0436-09 1.2196-06 2.9046-34 6.2176-03 4.509E-03 1.0426-06 

0.900006-01 5.009056-01 5.795696-01 1.125446 00 l .00049906 00 2.160176 00 1.005586 00 

9. 72TE-07 2.457E-02 1.094E-09 9.219E-04 2. 1976-04 6. 1946-03 4.565E-09 9.727E-07 

0.000006-01 4.740946-01 5.491186-01 1.175906 00 L. 19651096 00 2.160176 00 1.01574E 30 

0.0996-07 2.4596-02 1.9496-09 1.2196-04 2.0216-04 6.0496-01 4.5416-01 0.0996-07 

0.700006-01 4.449906-01 5.190096-01 1.22457E 00 1.19154426 00 2.L601TE 00 1.09045E 00 

0.1126-07 2.4416-02 2.0066-01 1.2196-04 1*8566-04 5. 9596-01 4.515E-09 8.1126-07 

0.600306-01 4.159496-01 4.915926-01 1.271456 00 1.24579446 00 2.160L7E 00 1.04940E 00 

7.909E-07 2.4696-02 2.066E-01 3.2196-04 1. 7006-04 5.0476-01 6.4876-03 7.30)6-07 

0.500306-01 9.07291E-0I 4.644506-01 1.116816 00 1.29927U6 00 2.140176 00 1.072456 00 

4.703E-07 2.4646-02 2.1296-09 9.2196-04 1.5546-04 5.7726-01 4. 4586-01 4.703E-07 

0.400006-01 1.404976-01 6.981726-01 1.140826 00 1.15291416 00 2.160176 00 1.099526 00 

6.0716-07 2.4686-02 2.1936-09 9.2196-04 1.4176-04 5.6716-01 4.6286-03 6.0716-07 


Figure 9 SAMPLE PRINTOUT FOR^EQU I LIBRIUM FLOW SOLUTION 
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DENSITY FIT-ALPHA- Z . t7 1 )6fr | E-0 l CONSTANT- 20699* JE-O? 


Figure 10 SAMPLE PRINTOUT OF DENSITY - FIT CONSTANTS FOR NONEQUILIBRIUM FLOW 
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NONEOUlLlBftIUM SOLUTION 


9. 900306-31 9.6)9206-31 9.56521E-01 3.38803E-01 *« 05961 85E-01 2.16017E 00 2.028216 00 -1.01630676 00 0 

0.990000 2.1 186-06 2.636E-0? 1.6506-03 3.2196-06 6.5I/E-U6 6.813E-0* 6.7566-01 2*1)86-0 6 

0.990000 2.1386-06 2.6)66-02 1.6S06-0) 3.2196-06 6.5776-06 6,8136-0) 6.756E-03 2.130E-O6 

9. 800306-01 8.906 1 36-01 9.106606-01 6.791176-01 6.35166 376-01 2.160176 00 1.503366 00 - 7 . 09667 766- 0 1 0 

0.980000 1.9926-06 2.6386-02 1. 6866-03 3.2196-06 6.2856-06 6.7756-03 6.7186-03 1.9926-06 

0.980000 1.9926-06 2. 6386-02 1.6866-01 3.2196-06 6.2856-06 6.7756-03 6.7386-0* I. 9926-06 

9.699996-01 8.393776-01 3.682996-01 5.867936-01 5.36017596-01 2.160176 00 1.287J6E 00 -5.36039786-01 0 

0.969999 1.8536-06 2.6636-02 1.5266-0* 3.2196-06 6.0066-06 6,7166-03 6.7236-03 1.8536-06 

0.970000 1.8536-06 2.6606-02 1.5266-01 3.2196-06 6.006E-06 6.7156-03 6.7236-0* 1.8536-06 

9.599996-01 7.907116-01 9.276516-01 6.775976-01 6.2263877E-01 2.160176 00 1.169866 00 -6. 121 38656-01 0 

0.959999 1.7216-06 2.6626-02 1.5666-01 3.2196-06 1.T39E-06 6.65*6-0* 6.7066-01 1.7216-06 

0.960000 1.7216-06 2.662E-02 1.5666-0* 3.2196-06 3.TJ9E-06 6.65*6-03 6.7066-03 1.7216-06 

9*699986-01 7.663206-01 7.88068E-01 7.57660E-01 7. 00362 70E - 01 2.L6017E 00 1.09858E 00 -3.13972596-01 0 

0*969998 1.5956-06 2.6666-02 1.6056-0* 3.2196-06 3.6866-06 6.5876-03 6.689E-03 1.5956-06 

0.950000 1.5956-06 2.6666-02 1.6056-03 3*2196-06 3.6866-06 6.587E-03 6.6896-03 1.593E-06 

9.399986-01 7.0011*6-01 7.50066E-01 8.30062E-01 7. 71982686-01 2.1631TE 00 1.053566 00 -2.31379606-01 0 

0.939998 1*6766-06 2.6666-02 1.6686-03 3.2196-06 3.2616-06 6.5196-03 6.6716-03 1.6766-06 

0.960000 1*6766-06 2.6666-02 1.6686-03 3.2196-06 3.2616-06 6.5196-01 6. 6716-01 1.676E-06 

9. 299976-01 6.580026-0L 7.131976-01 8.967296-01 8. 3916685E-0I 2.160176 00 1.025316 00 -i. 59087536-01 0 

0.929997 1.3666-06 2. 6696-02 1.6936-03 3.2196-06 3.010E-06 6.6686-03 6. 6526-03 1.3666-06 

0.930000 1 .1636-06 2.6696-02 1.6936-03 3.2196-06 3.009E-06 6. 6686-03 6. 6526-03 1.3616-06 

9. 199966-01 6.179066-01 6.780666-01 9.58B59E-01 9. 029*3296-01 2.160176 00 1.008866 00 -9.60096126-02 0 

0.919996 1.2576-06 2. 6516-02 I.T60E-0* 3.2196-06 2.7906-06 6.3766-03 6.6326-01 1.257E-06 

0.920000 1.2576-06 2.651E-02 1.7606-0* 3. 2196-06 2.790E-06 6.3766-03 6.6326-03 1.2576-06 

8. 96938E-01 5.257576-01 5.95056E-01 L.09965E 00 1.05260086 00 2.160176 00 1.002396 00 6.23310666-02 0 

0.896938 1.0176-06 2.6566-02 1.8676-03 3.2196-06 2.2886-06 6.1766-03 6.5776-03 1.0176-06 

0.896966 1.0166-06 2.6566-02 1.8676-0* 3.2196-06 2.288E-06 6.1766-03 6.5776-03 1.0L6E-06 

8.869366-01 6. 921076-01 5.66030E-01 1.151256 00 1.109129LE 00 2.160176 00 L.91012E DO 1.28376696-01 0 

0.886936 9.3086-07 2.6586-02 1.9216-03 3.2196-06 2.1076-06 6.092E-03 6.5536-03 9. 308E -07 

0.886946 9.2986-07 2.6586-02 1.9216-03 3.2196-06 2.1076-06 6.091E-O3 6.5536-03 9.2986-07 

8.76936E-0L 6.601376-01 5.361586-01 1.200776 00 1-1666765E 00 2.160176 00 1.022636 00 1.91938006-01 0 

0.876936 8.5036-07 2.6606-02 1.9786-03 3.2196-06 1.9366-06 6.0066-03 6. 5286-03 8.5036-07 

0.876946 8.6926-07 2.6406-02 1.9786-03 3.2196-06 1.9366-06 6.0066-03 6.5286-03 8.4926-07 

8.649326-01 4.297816-01 5.054056-01 1.268676 00 1.21928666 00 2.160176 00 1.039516 00 2.53639066-01 0 

0.866932 7.7516-07 2.6626-02 2.0366-03 3.2196-06 1.7766-06 5.9166-03 6.5016-03 7.7516-07 

0.866946 7.7*76-07 2.6626-02 2.0)66-0* 3.2196-06 1.7766-06 5.9136-03 6.5016-03 7.7376-07 

8.569286-01 6.009786-01 6.777626-01 1.29656E 00 1.273LS36E 00 2.16017E 00 1.060556 00 3. 13982656-01 0 

0.856928 7.0686-07 2.666E-02 2.097E-0I 3.2196-06 l. 6256-06 5.B20E-03 6.6736-03 T.068E-07 

0.B56946 7.0336-07 2.6666-02 2.0986-0* *.2196-06 1.625E-06 5.8196-0* 6.673E-03 T.033E-07 


Figure lla SAMPLE PRINTOUT FOR NONEQUILIBRIUM FLOW SOLUTION 
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*IQM£QUlLl§*lUiM SOLUTION 


8.4*925E-Ol 3.73667E-01 4.5ll*l£-0l t.))92?£ 00 l.J2645)'»E 00 2,16017* 00 1.0856)6 00 3*7)389726-01 0 

0.8**925 6.3956-07 2.467E-02 2.1616-0) 3. 2196-04 1.4046-04 5,7236-03 4,6636-03 6.395E-0T 

0*844444 6.378E-07 2.467E-02 2.1616-03 3.2196-04 l.4a3E-04 5.7226-03 4.44)6-0) 6.3T8E-07 

8. 349206-01 ). 477916-01 4.25575E-01 1.30Z61E 00 1.3793317E 00 2.160176 00 L.1L474E 00 *. 322 16 39E-0 1 0 

0.834920 5.7876-07 2.4696-0 2 2.2266-03 3.219E-04 1.3516-04 5.62)6-03 4.412E-03 5.707E-O7 

0.834944 5.7696-07 2.4696-02 2.2276-03 3.2196-04 1.35IE-04 5.6226-03 4.4UE-03 5.7696-07 

8.249146-01 3.23297E-01 4.0I018E-01 1.42484E 00 l. 43191376 00 2.160176 00 1.I4T94E 00 4.90777506-01 0 

0.824914 5.2256-07 2. 4716-02 2.2956-03 3.2196-04 1.2286-06 5.5196-03 4.1796-03 5.2256-07 

0.824944 5.2046-07 2.4716-02 2.2996-03 3.2196-0* 1.2286-0* 5*5196-03 4.378E-03 5.204E-07 

8.149076-01 3.001296-01 3.774476-01 1.466036 00 1.4943116E 00 2.160176 00 1.18536E 00 5.49)55076-01 0 

0.814907 4.7096-07 2.4746-02 2.3656-0) 3.2196-04 1.1136-0* 5.4136-0) 4.3446-03 4.705E-D7 

0.814944 4.6826-07 2.4746-02 2.3666-03 3.2196-04 1.11)6-04 5.4126-03 4.34)6-0) 4.682E-07 

8.048996-01 2.782)96-01 3.548396-01 1.50626E 00 1.5)662*06 OO 2.16017E 00 1.227206 00 6.08205)66-01 0 

0.804899 4.2256-07 2.476E-02 2.4)96-03 3.2196-04 1.306E-0* 5.30)6-03 4.307E-0) 4.2256-07 

0.804944 4.2006-07 2.4766-02 Z.440E-03 3.2196-04 1.006E-0* 5.3026-0) 4.3066-0) 4.200E-07 

7*948986-01 2.575766-01 3.33171E-01 l.5*56*E 00 1.5889399E 00 2.16017E 00 1.27371E 00 6.6757068E-0 l 0 

0.794988 3.7846-07 2.4786-02 2.9156-03 3.2196-0* 9.07*6-05 5.190E-0) *.26?E-03 3.784E-07 

0.794944 3,7576-07 2. 4786-02 2.5166-0) 3.2196-04 9.070F-05 5.1896-0) 4.266E-0) 3.757E-07 

7.948766-01 2.390956-01 3.12*2*6-01 t. 58*216 00 U6*l))88E 00 2.16Q17E 00 1.32523E 00 7.276812)6-01 0 

0.794876 3.3806-07 2.4916-02 2.59*6-0) 3.2196-0* 8.160E-05 5.07*6-0) *.2266-03 3.3806-07 

0.78*9** 3.3506-07 2. *816-02 2.5956-0) 3.2196-0* 8. 1556-05 5.072E-0) *.225E-0) 3.350E-0T 

7.748606-01 2.197486-01 2.92576E-01 1.6220SE 00 l. 69)89276 00 2.16017E 00 1.382116 00 7. 88T59566-01 0 

0.774860 3*0 L0E-0 7 2. *8)6-02 2.6T5E-0) 3.219E-04 7.3166-05 4.955E-0) 4.183E-0) 3.0106-07 

0.774944 2.978E-07 2. *8*6-02 2.6776-03 3-2196-0* 7.3126-05 *.95)6-0) *.101E-O) 2.978E-0T 

7.648276-01 1 • 9457)6-01 2.62*106-01 1.66572E 00 1.75001416 00 2.160176 00 1.50056E 00 9.0276632E-01 0 

0.764827 2.5886-07 2.4866-02 2.6736-Oi 3.2196-04 6.B00E-05 4.9696-0) 4.17)6-03 2.5006-07 

0.764944 2.5516-07 2.4046-02 2.6766-0) 3.2196-04 6.794E-05 4.9666-0) 4.1716-0) 2.5516-07 

7.548026-01 1.7908)6-01 2.451806-01 1.701986 00 1.8028355E 00 2.160176 00 1.571046 00 9.64907886-01 0 

0.754902 2.2926-07 2.4876-02 2,7606-03 3.219E-04 6.0576-05 4.9)96-0) 4.1286-0) 2.2926-07 

0.754944 2.25)6-07 2.487E-02 2.784E-0) 3.2196-04 6.0516-05 4.8)66-0) 4.1256-0) 2.2536-07 

7.447716-01 l *645446-01 2.297336-01 1. 737666 00 1.8560131E 00 2.160176 00 1.65026E 00 1.0289479E 00 0 

0.744771 2.024E-07 2.4906-02 2.8516-03 3.2196-04 5.3776-05 4.707E-0) 4.080E-0) 2.0246-07 

0.744944 1 .9826-07 2. 4906-02 2.855E-03 3.2196-04 5.3716-05 4. 7Q2E-0) 4.0766-0) 1.9026-07 

7.347)36-01 1.504206-01 2.l)059fc~01 L.7728LE 00 1.9095953E 00 2.163176 00 1.73654E 00 1.09507856 00 0 

0.7)473) 1.7826-0? 2.492E-02 2.9*36-0) 3.2196-0* *.7596-05 *.5716-03 4.0)06-0) 1.7026-07 

0.7)494* 1.7376-07 2.49)6-02 2.9406-0) 1.2196-04 *.7526-05 *.566E-0) *.0256-0) 1.7)76-07 

7.2*6876-01 1.3817)6-01 1.98l**fc-0l l.807*6E 00 1.96)62606 00 2.160176 00 L.83146E 00 1.163507*E 00 3 

0.724687 1.5646-07 2.49SE-02 3.039E-0) 3.2196-04 5.196E-05 4.*))E-0) 3.977E-03 1.56*6-07 

0.72*944 1.5176-07 Z. 4966-02 3.0456-03 3.2136-0* *.109E-O5 *.*276-0) 3.97U-0J 1.517E-07 
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I ih ii i ill ii li mi « in i mi mi 


Hill ilHIIIIIIliill 


miKjMiii 111 mi mi i mi mu 


hi i in ii idiii 


1111111111111 


i ilium is iismui i m i min i mi ihi 


N0HfcQUlLl8*lU« SOLUM (JM 


7.148296-01 
0. 714629 
0. 714944 

7.045*56-01 
0. 704565 
0. 70% 944 

6.944896*01 

0.694489 

0.494944 

6.84)946-01 

0.694)94 

Q.6B4944 

6. 742776-01 
0.674277 
0.474944 

6.641)27-01 

0 . 6641)2 

6.6)11)7 -01 
0. 66)1 15 

6 . 620996-0 1 
0.662099 

6.61 08)6-01 
0.66108) 

6.59971E-01 

0.6)9971 

6.598696-01 

0.6)8669 

6.577696-01 

0.6)7769 

6.566906-01 

0.6)6690 

6.5)48)6-01 

0.6)5485 

6.54104E-01 

0.654)04 


1 .262696 
1.3696-07 
1 • ) l 9fc -07 

1.151696- 

1.19)6-07 

1.1416-OT 

1.049)86- 

1.0) 96-07 
9.8)76-08 

9.524466 

9.0176-08 

8.4) 76-08 

8.6)5286 

7.8086-08 

7.2016-08 

7.805)06- 

6.7496-08 

7.725786 

6.6426-08 

7.647196 

6.5) 76-08 

7.569766 

6.4)66-08 

7.4857)6- 
6.)) 36-08 

7.40)146- 

6.2276-08 

7.321696- 
6.1 306-08 

7.241616- 

6.0) 06-09 

7.154816- 

5.9276-08 

7.0696)6- 

5.8226-08 


-01 


•02 1.4574 

2.5086-02 


•02 1 . )4)8 

2.5116-02 
2.512E-02 


■02 

2 . 

■02 

2. 

-02 


02 1.18)7 

2.517E-02 


2.5186-02 
02 1.1)9) 


2.018147)6 00 

2.160176 00 

04 3.6876-05 

4.29)6-0) 

04 3.6796-05 

4.2866-01 

2-07)19)46 00 

2.160176 00 

04 ). 2276-05 

4.1506-0) 

04 3.2186-05 

4.1416-03 

2.12879796 00 

2.160176 00 

04 2.8146-05 

4.004E-0) 

■04 2.8046-05 

1.9946-0) 

2.18498986 00 

2.160176 00 

■04 2.44)6-05 

3.B57E-0) 

-04 2.4)36-05 

3.845E-03 

2.24179446 00 

2.160176 00 

-04 2.1126-05 

). 7086-0) 

-04 2.1016-05 

3.69)6-0) 

2.30018686 00 

2.160246 00 

-04 1.8186-05 

*.5586-0) 

2. )05 7093E 00 

?• 16025b 00 

-04 l . 791 E -05 

3.543E-0) 

2.30989166 00 

2.160246 00 

-04 1.76)6-05 

). 5286-0) 

2.31571626 00 

2.160256 00 

-04 1.7)66-05 

3.5136-0) 

2.J216189E 00 

2*160256 00 

-04 1.7076-05 

3.496E-0) 

2.32824216 00 

2.160256 00 

-04 1.679E-05 

3.480E-0) 

2.3)4049)6 00 

2.160256 00 

-04 l. 6506-05 

). 46)6-0) 

2.34053016 00 

2.160256 00 

-04 1.62)6-05 

3.447E-0) 

2.34700776 00 

2.160256 00 

-04 1.5946-05 

1.4296-03 

2.35419796 00 

2.160256 00 

-04 1.5656-05 

3.4126-0) 


1.9)59)6 00 1.2)445626 00 

9.9226-0) 1.3696-07 

3.9146-0) l. *196-07 

2.051066 00 1. 90871406 00 

3.8646-0) 1.1956-07 

1.8)46-0) 1.141E-07 

2.17796E 00 1. >877)896 00 

3.8046-0) 1.0)96 -07 

3.7916-0) 9*8)76-08 

2.315006 00 1.471882)6 00 

3.741E-0) 9. 0176-08 

3.726E-0) 8.4)76-08 

2.472706 00 1.56148)36 00 

3.6766-0) 7.8086-08 

3.6576-0) 7 .2016-08 

2.64)796 00 1.65691216 00 

3.609E-0) 6. T496-08 

2.662096 00 1.6669L216 00 

3.602E-0) 6.642E-08 

2.68047E 00 1.6769121E 00 

3.5)56-0) 6.5)76-08 

2.698926 00 1.68691216 00 

). 5886-0) 6.4)66-08 

2.719)06 00 1.59T9121E 00 

3.580E-0) 6.1)36-08 

2.7)9766 00 1.70891206 00 

3.5726-0) 6.2276-08 

2. 760) IE 00 1.71991206 00 

3.5646-0) 6.1)06-08 

2.T8095E 00 1.7)091206 00 

3.5576-0) 6.0)06-08 

2.80375E 00 1.74)01206 00 

3.5496-0) 5.9276-08 

2. 826656 00 I.TS51120E 00 

3.540E-0) 5.8226-08 


0 


0 


0 


0 


0 


1 


1 

l 


l 


L 

1 


l 


1 

1 

l 


Figure lie SAMPLE PRINTOUT FOR NONEQUILIBRIUM FLOW SOLUTION 



NASA- Langley, 1966 


4 ^ 

tv 


1000 0. 714944 


-0* 

5.26-06 

-1.56-02 

-1.56-02 

2.26-07 

1.16-02 

-l. 16-02 

-1. 16-02 

-0. 

7.16-01 

-9.16-01 

-9.16-01 

4. 06-07 

1 .46-02 

-8.66-01 

-8.56-03 

1000 0. 

704944 



-2. 26-37 

J.5E-06 

-1.96-02 

-1.96-02 

2.26-37 

8.96-01 

-1.46-02 

-1.56-02 

-0. 

5. 16-03 

-1*16-02 

-1.16-02 

-0. 

l. 16-02 

-l. 16-02 

-l. 16-02 

1000 3.< 

694944 



-0. 

2.1E-06 

-2.46-02 

-2.16-02 

2.26-07 

7.1E-0J 

-1.76-02 

-1.76-02 

-0. 

1. 76-01 

-1.46-02 

-1.46-02 

-0* 

3.96-01 

-1.16-02 

-l. 16-02 

1000 O.i 

584944 



-o.. 

1.56-06 

-2. 96-02 

-2.96-02 

-0. 

6.06-03 

-2.16-02 

-2. 16-02 

-0. 

2 .66-03 

-1.76-02 

-1.76-02 

-0. 

7.06-03 

-1.66-02 

-L . 56-02 

1000 0.674944 



2.26-07 

9.16-07 

-1.76-02 

-1. 46-02 

-0. 

4.9E-01 

-2.66-02 

-2.86-02 

-0. 

1.96-01 

-2.16-02 

-2.16-02 

4.86-0. 

5.SE-01 

-2.16-02 

-2.06-02 

1000 0. 6641 >2 



-4.76-02 

5.76-0? 

-2.66-08 

-l. 16-07 

-1.16-32 

1. 86-01 

-1.26-06 

0. 

-2.66-02 

l. >6-01 

-1.46-05 

-7.56-04 

-2.66-02 

4.26-01 

-1.16-04 

0. 

1000 0.66*115 



-4.8E-02 

5.46-07 

-2.66-05 

-1. 16-07 

-1. 06-02 

1.76-01 

-1.2E-04 

0. 

-2.76-02 

1.26-01 

-1. 16-05 

-7.26-04 

-2.66-02 

4. 16-03 

-I.IE-04 

0. 

1000 0.662009 



-4.86-02 

5.26-07 

-2.56-08 

-1.06-07 

-1.46-02 

1.6E-03 

-l. 26-04 

0. 

-2.76-02 

1.26-01 

-1. 26-05 

-7.26-04 

-2.76-02 

1.9E-01 

-1.1E-04 

0. 

1000 0.661035 



-4.96-02 

4.96-07 

-2.46-08 

-9.96-08 

-1.46-02 

1.5E-01 

-1.26-04 

0. 

-2.8E-02 

1.26-01 

-1.26-05 

-7.16-04 

-2.76-02 

1. BE -01 

-1.06-04 

0. 

1000 0.659971 



-5.06-02 

4.76-07 

-2.16-08 

-9. 76-08 

-1.56-02 

J.5E-01 

-1.26-04 

0. 

-2 .86-02 

l. 16-01 

-1. 16-05 

-7.06-04 

-2.86-02 

1.76-01 

-1.06-04 

0. 

1000 0.658869 



-5.16-32 

4 .46-07 

-2.16-09 

-9.46-08 

-1.66-02 

1. 46-03 

-1.26-04 

0 . 

-2.96-02 

L. IE-01 

- J. 16-05 

-6.9fc 05 

-2.86-02 

1. 66-03 

-l. 06-05 

0 . 

1000 0.657769 



-5.26-02 

4.2E-J7 

-2.26-09 

-9. 16-08 

-1.76-02 

1. 16-01 

-1.2E-04 

0. 

-1.06-02 

t .06-03 

-i. 16-05 

"8. Jt-U4 

-2.96-02 

1.56-03 

-l. 06-04 

0. 

1000 0.656680 



-5.46-02 

4.0E-07 

-2.16-08 

-9.06-08 

-1.8E-02 

1.26-01 

-1.76-05 

0 . 


-o. 

2.26-07 

4.86-07 

-0. 

1.96-08 

4.06-02 

1.56-01 

1.46-08 

-1.5E-02 
-1. 16-02 
-2.36-03 
7.16-02 

-1.56-02 

-1.16-02 

-2.16-03 

7.76-02 

-2.2E-07 

2.26-07 

4.86-07 

2.26-07 

2.76-08 
1.06-07 
L. 16-01 
Z.3E-08 

-1.96-02 

-1.46-02 

-2.86-01 

8.76-02 

-1.96-02 

-1.46-02 

-2.86-03 

9.36-02 

-0. 

2.26-07 

2.26-07 

-0. 

1. 16- 08 

2.16- 02 
8.96-02 
1.66-06 

-2.5E-02 

-1.76-02 

-1.46-01 

1.06-01 

-2. 36-02 
-1.76-02 
-1.46-03 
1 .16-01 

- 3. 

-0. 

-0. 

2.26-07 
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Figure 12 SAMPLE PRINTOUT OF REACTION RATES COMPUTED FOR NONEQUILIBRIUM FLOW 



